New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 2905 for branches/2011/dev_r2739_LOCEAN8_ZTC – NEMO

Ignore:
Timestamp:
2011-10-12T14:28:01+02:00 (13 years ago)
Author:
mlelod
Message:

first commit, compilation ok, see ticket/863?

Location:
branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM
Files:
26 edited

Legend:

Unmodified
Added
Removed

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/GYRE/EXP00/namelist

    r2715 r2905  
    500500!!                      ***  Dynamics namelists  *** 
    501501!!====================================================================== 
     502!!   nam_vvl       vertical coordinate options 
    502503!!   namdyn_adv    formulation of the momentum advection 
    503504!!   namdyn_vor    advection scheme 
    … …  
    507508!!====================================================================== 
    508509! 
     510!----------------------------------------------------------------------- 
     511&nam_vvl    !   vertical coordinate options 
     512!----------------------------------------------------------------------- 
     513&nam_vvl 
     514   ln_vvl_zstar  = .false. !  zstar vertical coordinate                     
     515   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations 
     516   ln_vvl_layer  = .false. !  full layer vertical coordinate 
     517   ahe3          = 0.e0    !  thickness diffusion coefficient  
     518/ 
    509519!----------------------------------------------------------------------- 
    510520&namdyn_adv    !   formulation of the momentum advection 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist

    r2735 r2905  
    500500!!                      ***  Dynamics namelists  *** 
    501501!!====================================================================== 
     502!!   nam_vvl       vertical coordinate options 
    502503!!   namdyn_adv    formulation of the momentum advection 
    503504!!   namdyn_vor    advection scheme 
    … …  
    507508!!====================================================================== 
    508509! 
     510!----------------------------------------------------------------------- 
     511&nam_vvl    !   vertical coordinate options 
     512!----------------------------------------------------------------------- 
     513&nam_vvl 
     514   ln_vvl_zstar  = .false. !  zstar vertical coordinate                     
     515   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations 
     516   ln_vvl_layer  = .false. !  full layer vertical coordinate 
     517   ahe3          = 0.e0    !  thickness diffusion coefficient  
     518/ 
    509519!----------------------------------------------------------------------- 
    510520&namdyn_adv    !   formulation of the momentum advection 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml

    r2729 r2905  
    4848   <field id="botpres"      description="Pressure at sea floor"                     unit="dbar"                     /> 
    4949   <field id="cellthc"      description="Cell thickness"                            unit="m"     axis_ref="deptht"  /> 
     50   <!-- variables available with key_vvl --> 
     51   <field id="e3t_n"        description="vertical scale factor"                     unit="m"                        /> 
     52   <field id="e3tdef"       description="vertical scale factor variance"            unit="%^2"                      /> 
     53   <field id="dept_n"       description="depth"                                     unit="m"                        /> 
    5054     </group> 
    5155 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_ar5.xml

    r2561 r2905  
    4848   <field id="botpres"      description="Pressure at sea floor"                     unit="dbar"                     /> 
    4949   <field id="cellthc"      description="Cell thickness"                            unit="m"     axis_ref="deptht"  /> 
     50   <!-- variables available with key_vvl --> 
     51   <field id="e3t_n"        description="vertical scale factor"                     unit="m"                        /> 
     52   <field id="e3tdef"       description="vertical scale factor variance"            unit="%^2"                      /> 
     53   <field id="dept_n"       description="depth"                                     unit="m"                        /> 
    5054     </group> 
    5155 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

    r2715 r2905  
    500500!!                      ***  Dynamics namelists  *** 
    501501!!====================================================================== 
     502!!   nam_vvl       vertical coordinate options 
    502503!!   namdyn_adv    formulation of the momentum advection 
    503504!!   namdyn_vor    advection scheme 
    … …  
    507508!!====================================================================== 
    508509! 
     510!----------------------------------------------------------------------- 
     511&nam_vvl    !   vertical coordinate options 
     512!----------------------------------------------------------------------- 
     513&nam_vvl 
     514   ln_vvl_zstar  = .true.  !  zstar vertical coordinate                     
     515   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations 
     516   ln_vvl_layer  = .false. !  full layer vertical coordinate 
     517   ahe3          = 0.e0    !  thickness diffusion coefficient  
     518/ 
    509519!----------------------------------------------------------------------- 
    510520&namdyn_adv    !   formulation of the momentum advection 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm

    r2670 r2905  
    1  bld::tool::fppkeys key_trabbl key_vectopt_loop key_orca_r2 key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_dtatem key_dtasal key_tradmp key_zdftke key_zdfddm key_zdftmx key_iomput  
     1 bld::tool::fppkeys key_trabbl key_vectopt_loop key_orca_r2 key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_dtatem key_dtasal key_tradmp key_zdftke key_zdfddm key_zdftmx key_iomput key_vvl 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist

    r2715 r2905  
    500500!!                      ***  Dynamics namelists  *** 
    501501!!====================================================================== 
     502!!   nam_vvl       vertical coordinate options 
    502503!!   namdyn_adv    formulation of the momentum advection 
    503504!!   namdyn_vor    advection scheme 
    … …  
    508509!!====================================================================== 
    509510! 
     511!----------------------------------------------------------------------- 
     512&nam_vvl    !   vertical coordinate options 
     513!----------------------------------------------------------------------- 
     514&nam_vvl 
     515   ln_vvl_zstar  = .false. !  zstar vertical coordinate                     
     516   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations 
     517   ln_vvl_layer  = .false. !  full layer vertical coordinate 
     518   ahe3          = 0.e0    !  thickness diffusion coefficient  
     519/ 
    510520!----------------------------------------------------------------------- 
    511521&namdyn_adv    !   formulation of the momentum advection 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/CONFIG/POMME/EXP00/namelist

    r2650 r2905  
    500500!!                      ***  Dynamics namelists  *** 
    501501!!====================================================================== 
     502!!   nam_vvl       vertical coordinate options 
    502503!!   namdyn_adv    formulation of the momentum advection 
    503504!!   namdyn_vor    advection scheme 
    … …  
    507508!!====================================================================== 
    508509! 
     510!----------------------------------------------------------------------- 
     511&nam_vvl    !   vertical coordinate options 
     512!----------------------------------------------------------------------- 
     513&nam_vvl 
     514   ln_vvl_zstar  = .false. !  zstar vertical coordinate                     
     515   ln_vvl_ztilde = .false. !  hybrid verticalcoordinate: only high frequency variations 
     516   ln_vvl_layer  = .false. !  full layer vertical coordinate 
     517   ahe3          = 0.e0    !  thickness diffusion coefficient  
     518/ 
    509519!----------------------------------------------------------------------- 
    510520&namdyn_adv    !   formulation of the momentum advection 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90

    r2715 r2905  
    249249             
    250250            !  energy needed to bring ocean surface layer until its freezing 
    251             qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1)   & 
     251            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj)   & 
    252252                &          * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda ) 
    253253             

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r2715 r2905  
    139139   !! All coordinates 
    140140   !! --------------- 
    141    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_1           !: depth of T-points (sum of e3w) (m) 
    142    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_1, gdepw_1   !: analytical depth at T-W  points (m) 
    143    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v_1  , e3f_1     !: analytical vertical scale factors at  V--F 
    144    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_1  , e3u_1     !:                                       T--U  points (m) 
    145    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3vw_1             !: analytical vertical scale factors at  VW-- 
    146    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_1  , e3uw_1    !:                                       W--UW  points (m) 
    147    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_b              !: before         -      -      -    -   T      points (m) 
    148    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_b  , e3v_b     !:   -            -      -      -    -   U--V   points (m) 
     141   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_n           !: now depth of T-points (sum of e3w) (m) 
     142   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_n, gdepw_n   !: now depth at T-W  points (m) 
     143   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_n              !: now    vertical scale factors at  T       point  (m) 
     144   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_n  , e3v_n     !:            -      -      -    -   U --V   points (m) 
     145   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_n  , e3f_n     !:            -      -      -    -   w --F   points  (m) 
     146   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3uw_n , e3vw_n    !:            -      -      -    -   UW--VW  points (m) 
     147   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_b              !: before     -      -      -    -   T       points (m) 
     148   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_b  , e3v_b     !:   -        -      -      -    -   U --V   points (m) 
     149   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3uw_b , e3vw_b    !:   -        -      -      -    -   UW--VW  points (m) 
     150   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_a              !: after      -      -      -    -   T       point  (m) 
     151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3u_a  , e3v_a     !:   -        -      -      -    -   U --V   points (m) 
    149152#else 
    150153   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag 
    151154#endif 
    152    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur  , hvr    !: inverse of u and v-points ocean depth (1/m) 
    153    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu   , hv     !: depth at u- and v-points (meters) 
    154    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0 , hv_0   !: refernce depth at u- and v-points (meters) 
     155   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur    , hvr       !: inverse of u and v-points ocean depth (1/m) 
     156   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu     , hv        !: depth at u- and v-points (meters) 
     157   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0               !: refernce depth at t-       points (meters) 
     158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0   , hv_0      !: refernce depth at u- and v-points (meters) 
     159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1ur   , e2vr      !: scale factor coeffs at U--V points 
     160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12t   , e12t_1    !: horizontal cell surface and inverse  at T    points 
     161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12u   , e12u_1    !: horizontal cell surface and inverse  at U    points 
     162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e12v_1 , e12f_1    !: inverse horizontal cell surface      at V--F   points 
    155163 
    156164   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
    … …  
    277285         ! 
    278286#if defined key_vvl 
    279       ALLOCATE( gdep3w_1(jpi,jpj,jpk) , e3v_1(jpi,jpj,jpk) , e3f_1 (jpi,jpj,jpk) ,                           & 
    280          &      gdept_1 (jpi,jpj,jpk) , e3t_1(jpi,jpj,jpk) , e3u_1 (jpi,jpj,jpk) ,                           & 
    281          &      gdepw_1 (jpi,jpj,jpk) , e3w_1(jpi,jpj,jpk) , e3vw_1(jpi,jpj,jpk) , e3uw_1(jpi,jpj,jpk) ,     & 
    282          &      e3t_b   (jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk)                       , STAT=ierr(5) ) 
    283 #endif 
    284          ! 
    285       ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) ,     & 
    286          &      hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , STAT=ierr(6) ) 
     287      ALLOCATE( gdep3w_n(jpi,jpj,jpk) , e3t_n (jpi,jpj,jpk) , e3u_n (jpi,jpj,jpk) ,                           & 
     288         &      gdept_n (jpi,jpj,jpk) , e3v_n (jpi,jpj,jpk) , e3w_n (jpi,jpj,jpk) ,                           & 
     289         &      gdepw_n (jpi,jpj,jpk) , e3f_n (jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , e3uw_n(jpi,jpj,jpk) ,     & 
     290         &      e3t_b   (jpi,jpj,jpk) , e3u_b (jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) ,                           & 
     291         &      e3uw_b  (jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,                                                 & 
     292         &      e3t_a   (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk)                       , STAT=ierr(5) ) 
     293#endif 
     294         ! 
     295      ALLOCATE( hu  (jpi,jpj) , hur   (jpi,jpj) , hu_0  (jpi,jpj) , ht_0  (jpi,jpj) ,     & 
     296         &      hv  (jpi,jpj) , hvr   (jpi,jpj) , hv_0  (jpi,jpj) ,                       & 
     297         &      e1ur(jpi,jpj) , e2vr  (jpi,jpj) , e12t  (jpi,jpj) , e12t_1(jpi,jpj) ,     & 
     298         &      e12u(jpi,jpj) , e12u_1(jpi,jpj) , e12v_1(jpi,jpj) , e12f_1(jpi,jpj) , STAT=ierr(6)  ) 
    287299         ! 
    288300      ALLOCATE( gdept_0(jpk) , gdepw_0(jpk) ,                                     & 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r2528 r2905  
     1 
    12MODULE domain 
    23   !!============================================================================== 
    … …  
    8182                             CALL dom_zgr      ! Vertical mesh and bathymetry 
    8283                             CALL dom_msk      ! Masks 
    83       IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh 
     84      IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh 
    8485      ! 
    8586      IF( lk_c1d ) THEN                        ! 1D configuration  
    … …  
    8990      END IF 
    9091      ! 
     92      ! - ML - only used in domvvl but could be usefull in many routines 
     93      e12t  (:,:) = e1t(:,:) * e2t(:,:) 
     94      e12u  (:,:) = e1u(:,:) * e2u(:,:) 
     95      e12u_1(:,:) = 0.5 / e12u(:,:) 
     96      e12v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) ) 
     97      e12f_1(:,:) = 0.5 / ( e1f(:,:) * e2f(:,:) ) 
     98      ! - ML - used in domvvl and traldf_(lab/bilap/iso) 
     99      e1ur  (:,:) = e2u(:,:) / e1u(:,:) 
     100      e2vr  (:,:) = e1v(:,:) / e2v(:,:) 
     101      ! 
    91102      hu(:,:) = 0.e0                           ! Ocean depth at U- and V-points 
    92103      hv(:,:) = 0.e0 
    93104      DO jk = 1, jpk 
    94          hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    95          hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
     105         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 
     106         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
    96107      END DO 
    97108      !                                        ! Inverse of the local depth 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r2715 r2905  
    66   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code 
    77   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate 
    8    !!---------------------------------------------------------------------- 
     8   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: 
     9   !!                                          vvl option includes z_star and z_tilde coordinates 
    910#if defined key_vvl 
    1011   !!---------------------------------------------------------------------- 
    1112   !!   'key_vvl'                              variable volume 
    1213   !!---------------------------------------------------------------------- 
    13    !!   dom_vvl     : defined coefficients to distribute ssh on each layers 
    1414   !!---------------------------------------------------------------------- 
     15   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
     16   !!   dom_vvl_nxt      : Compute next vertical scale factors 
     17   !!   dom_vvl_swp      : Swap vertical scale factors and update the vertical grid 
     18   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 
     19   !!   dom_vvl_rst      : read/write restart file 
     20   !!   dom_vvl_ctl      : Check the vvl options 
     21   !!---------------------------------------------------------------------- 
     22   !! * Modules used 
    1523   USE oce             ! ocean dynamics and tracers 
    1624   USE dom_oce         ! ocean space and time domain 
    17    USE sbc_oce         ! surface boundary condition: ocean 
    18    USE phycst          ! physical constants 
     25   USE sbc_oce         ! ocean surface boundary condition 
    1926   USE in_out_manager  ! I/O manager 
     27   USE iom             ! I/O manager library 
     28   USE restart         ! ocean restart 
    2029   USE lib_mpp         ! distributed memory computing library 
    2130   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    2433   PRIVATE 
    2534 
    26    PUBLIC   dom_vvl       ! called by domain.F90 
    27    PUBLIC   dom_vvl_alloc ! called by nemogcm.F90 
    28  
    29    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ee_t, ee_u, ee_v, ee_f   !: ??? 
    30    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: ???  
    31  
    32    REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra  
    33       !                                                              ! except at nit000 (=rdttra) if neuler=0 
     35   !! * Routine accessibility 
     36   PUBLIC dom_vvl_init       ! called by domain.F90 
     37   PUBLIC dom_vvl_sf_nxt     ! called by step.F90 
     38   PUBLIC dom_vvl_sf_swp     ! called by step.F90 
     39   PUBLIC dom_vvl_interpol   ! called by dynnxt.F90 
     40 
     41   !!* Namelist nam_vvl 
     42   LOGICAL , PUBLIC                                      :: ln_vvl_zstar  = .TRUE.    ! zstar  vertical coordinate 
     43   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde = .FALSE.   ! ztilde vertical coordinate 
     44   LOGICAL , PUBLIC                                      :: ln_vvl_layer  = .FALSE.   ! level  vertical coordinate 
     45   LOGICAL , PUBLIC                                      :: ln_vvl_kepe   = .FALSE.   ! kinetic/potential energy transfer 
     46   !                                                                                  ! conservation: not used yet 
     47   REAL(wp), PUBLIC                                      :: ahe3          =  0.e0     ! thickness diffusion coefficient 
     48 
     49   !! * Module variables 
     50   INTEGER                                               :: nvvl                      ! choice of vertical coordinate 
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td              ! thickness diffusion transport 
     52   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                   ! low frequency part of hz divergence 
     53   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_t_b, e3t_t_n, e3t_t_a ! baroclinic scale factors 
     54   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_restore_e3t           ! retoring period for scale factors 
     55   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_restore_hdv           ! retoring period for low freq. divergence 
     56 
    3457 
    3558   !! * Substitutions 
    … …  
    3760#  include "vectopt_loop_substitute.h90" 
    3861   !!---------------------------------------------------------------------- 
    39    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     62   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)  
    4063   !! $Id$ 
    4164   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4265   !!---------------------------------------------------------------------- 
    43 CONTAINS        
     66 
     67CONTAINS 
    4468 
    4569   INTEGER FUNCTION dom_vvl_alloc() 
    4670      !!---------------------------------------------------------------------- 
    47       !!                ***  ROUTINE dom_vvl_alloc  *** 
    48       !!---------------------------------------------------------------------- 
    49       ! 
    50       ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) ,     & 
    51          &      ee_t(jpi,jpj)     , ee_u(jpi,jpj)     , ee_v(jpi,jpj)     , ee_f(jpi,jpj)     ,     & 
    52          &      r2dt        (jpk)                                                             , STAT=dom_vvl_alloc ) 
    53          ! 
    54       IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
    55       IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 
    56       ! 
     71      !!                ***  FUNCTION dom_vvl_alloc  *** 
     72      !!---------------------------------------------------------------------- 
     73      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     74         ALLOCATE(                                                                                             & 
     75            &      un_td  (jpi,jpj,jpk) , vn_td  (jpi,jpj,jpk) , hdiv_lf(jpi,jpj,jpk) ,                        & 
     76            &      e3t_t_b(jpi,jpj,jpk) , e3t_t_n(jpi,jpj,jpk) , e3t_t_a(jpi,jpj,jpk) ,                        & 
     77            &      frq_restore_e3t(jpi,jpj) , frq_restore_hdv(jpi,jpj)                , STAT= dom_vvl_alloc ) 
     78         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
     79         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 
     80      ENDIF 
     81 
    5782   END FUNCTION dom_vvl_alloc 
    5883 
    5984 
    60    SUBROUTINE dom_vvl 
    61       !!---------------------------------------------------------------------- 
    62       !!                ***  ROUTINE dom_vvl  *** 
     85   SUBROUTINE dom_vvl_init 
     86      !!---------------------------------------------------------------------- 
     87      !!                ***  ROUTINE dom_vvl_init  *** 
    6388      !!                    
    64       !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread 
    65       !!               ssh over the whole water column (scale factors) 
    66       !!---------------------------------------------------------------------- 
    67       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    68       USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3     ! 2D workspace 
    69       ! 
    70       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    71       REAL(wp) ::   zcoefu , zcoefv   , zcoeff                   ! local scalars 
    72       REAL(wp) ::   zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1   !   -      - 
    73       !!---------------------------------------------------------------------- 
    74  
    75       IF( wrk_in_use(2, 1,2,3) ) THEN 
    76          CALL ctl_stop('dom_vvl: requested workspace arrays unavailable')   ;   RETURN 
    77       ENDIF 
    78  
    79       IF(lwp) THEN 
    80          WRITE(numout,*) 
    81          WRITE(numout,*) 'dom_vvl : Variable volume initialization' 
    82          WRITE(numout,*) '~~~~~~~~  compute coef. used to spread ssh over each layers' 
    83       ENDIF 
    84        
    85       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl : unable to allocate arrays' ) 
    86  
    87       fsdept(:,:,:) = gdept (:,:,:) 
    88       fsdepw(:,:,:) = gdepw (:,:,:) 
    89       fsde3w(:,:,:) = gdep3w(:,:,:) 
    90       fse3t (:,:,:) = e3t   (:,:,:) 
    91       fse3u (:,:,:) = e3u   (:,:,:) 
    92       fse3v (:,:,:) = e3v   (:,:,:) 
    93       fse3f (:,:,:) = e3f   (:,:,:) 
    94       fse3w (:,:,:) = e3w   (:,:,:) 
    95       fse3uw(:,:,:) = e3uw  (:,:,:) 
    96       fse3vw(:,:,:) = e3vw  (:,:,:) 
    97  
    98       !                                 !==  mu computation  ==! 
    99       ee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level 
    100       ee_u(:,:) = fse3u_0(:,:,1) 
    101       ee_v(:,:) = fse3v_0(:,:,1) 
    102       ee_f(:,:) = fse3f_0(:,:,1) 
    103       DO jk = 2, jpkm1                          ! Sum of the masked vertical scale factors 
    104          ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 
    105          ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
    106          ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
    107          DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask 
    108             ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk) 
    109          END DO 
    110       END DO   
    111       !                                         ! Compute and mask the inverse of the local depth at T, U, V and F points 
    112       ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1) 
    113       ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1) 
    114       ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1) 
    115       DO jj = 1, jpjm1                               ! f-point case fmask cannot be used  
    116          ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 
     89      !! ** Purpose :  Initialization of all scale factors, depths 
     90      !!               and water column heights 
     91      !! 
     92      !! ** Method  :  - use restart file and/or initialize 
     93      !!               - interpolate scale factors 
     94      !! 
     95      !! ** Action  : - fse3t_(n/b) and e3t_t_(n/b) 
     96      !!              - Regrid: fse3(u/v)_n 
     97      !!                        fse3(u/v)_b        
     98      !!                        fse3w_n            
     99      !!                        fse3(u/v)w_b       
     100      !!                        fse3(u/v)w_n       
     101      !!                        fsdept_n, fsdepw_n and fsde3w_n 
     102      !!                        h(u/v) and h(u/v)r 
     103      !!              - ht_0 and ht1_0 
     104      !! 
     105      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
     106      !!---------------------------------------------------------------------- 
     107      !! * Local declarations 
     108      INTEGER ::   jk 
     109      !!---------------------------------------------------------------------- 
     110 
     111      IF(lwp) WRITE(numout,*) 
     112      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
     113      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     114 
     115#if defined key_zco 
     116      CALL ctl_stop( 'dom_vvl_init : options key_zco/key_dynspg_rl are incompatible with variable volume option key_vvl') 
     117#endif 
     118 
     119      ! choose vertical coordinate (z_star, z_tilde or layer) 
     120      ! ========================== 
     121      CALL dom_vvl_ctl 
     122 
     123      ! Allocate module arrays 
     124      ! ====================== 
     125      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
     126 
     127      ! read or initialize e3t_t_(b/n) and fse3t_(b/n)  
     128      ! ============================================== 
     129      CALL dom_vvl_rst( nit000, 'READ' ) 
     130 
     131      ! Reconstruction of all vertical scale factors at now and before time steps 
     132      ! ========================================================================= 
     133      ! Horizontal scale factor interpolations 
     134      ! -------------------------------------- 
     135      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 
     136      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 
     137      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 
     138      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
     139      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 
     140      ! Vertical scale factor interpolations 
     141      ! ------------------------------------ 
     142      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  ) 
     143      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 
     144      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 
     145      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 
     146      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 
     147      ! t- and w- points depth 
     148      ! ---------------------- 
     149      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1) 
     150      fsdepw_n(:,:,1) = 0.e0 
     151      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 
     152      DO jk = 2, jpk 
     153         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 
     154         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 
     155         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:) 
    117156      END DO 
    118       CALL lbc_lnk( ee_f, 'F', 1. )                  ! lateral boundary condition on ee_f 
    119       ! 
    120       DO jk = 1, jpk                            ! mu coefficients 
    121          mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels 
    122          muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk)     ! U-point at T levels 
    123          muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels 
    124       END DO 
    125       DO jk = 1, jpk                                 ! F-point : fmask=shlat at coasts, use the product of umask 
    126          DO jj = 1, jpjm1 
    127                muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels 
    128          END DO 
    129          muf(:,jpj,jk) = 0.e0 
    130       END DO 
    131       CALL lbc_lnk( muf, 'F', 1. )                   ! lateral boundary condition 
    132  
    133  
    134       hu_0(:,:) = 0.e0                          ! Reference ocean depth at U- and V-points 
     157      ! Reference water column height at t-, u- and v- point 
     158      ! ---------------------------------------------------- 
     159      ht_0(:,:) = 0.e0 
     160      hu_0(:,:) = 0.e0 
    135161      hv_0(:,:) = 0.e0 
    136162      DO jk = 1, jpk 
     163         ht_0(:,:) = ht_0(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 
    137164         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
    138165         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
    139166      END DO 
    140        
    141       ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 
    142       ! for ssh and scale factors 
    143       zs_t  (:,:) =       e1t(:,:) * e2t(:,:) 
    144       zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:) 
    145       zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:) 
    146  
    147       DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points 
    148          DO ji = 1, jpim1   ! NO vector opt. 
    149             zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj) 
    150             zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj) 
    151             zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) ) 
    152             ! before fields 
    153             zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  ) 
    154             zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  ) 
    155             zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1) 
    156             sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 
    157             sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 
    158             ! now fields 
    159             zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  ) 
    160             zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  ) 
    161             zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1) 
    162             zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1) 
    163             sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j ) 
    164             sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 ) 
    165             sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 ) 
    166          END DO 
     167 
     168   END SUBROUTINE dom_vvl_init 
     169 
     170 
     171   SUBROUTINE dom_vvl_sf_nxt( kt )  
     172      !!---------------------------------------------------------------------- 
     173      !!                ***  ROUTINE dom_vvl_sf_nxt  *** 
     174      !!                    
     175      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt, 
     176      !!                 tranxt and dynspg routines 
     177      !! 
     178      !! ** Method  :  - z_tilde_case: after scale factor increment computed with 
     179      !!                               high frequency part of horizontal divergence + retsoring to 
     180      !!                               towards the background grid + thickness difusion. 
     181      !!                               Then repartition of ssh INCREMENT proportionnaly 
     182      !!                               to the "baroclinic" level thickness. 
     183      !!               - z_star case:  Repartition of ssh proportionnaly to the level thickness. 
     184      !! 
     185      !! ** Action  :  - hdiv_lf: restoring towards full baroclinic divergence in z_tilde case 
     186      !!               - e3t_t_a: after increment of vertical scale factor  
     187      !!                          in z_tilde case 
     188      !!               - fse3(t/u/v)_a 
     189      !! 
     190      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
     191      !!---------------------------------------------------------------------- 
     192      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released 
     193      USE oce     , ONLY:   ze3t    => ta         ! ua used as workspace 
     194      USE wrk_nemo, ONLY:   zht     => wrk_2d_1   ! 2D REAL workspace 
     195      USE wrk_nemo, ONLY:   z_scale => wrk_2d_2   ! 2D REAL workspace 
     196      USE wrk_nemo, ONLY:   zwu     => wrk_2d_3   ! 2D REAL workspace 
     197      USE wrk_nemo, ONLY:   zwv     => wrk_2d_4   ! 2D REAL workspace 
     198      USE wrk_nemo, ONLY:   zhdiv   => wrk_2d_5   ! 2D REAL workspace 
     199      !! * Arguments 
     200      INTEGER, INTENT( in )                  :: kt                    ! time step 
     201      !! * Local declarations 
     202      INTEGER                                :: ji, jj, jk            ! dummy loop indices 
     203      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers 
     204      REAL(wp)                               :: z2dt                  ! temporary scalars 
     205      REAL(wp)                               :: z_def_max             ! temporary scalar 
     206      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars 
     207      LOGICAL                                :: ln_debug = .FALSE.    ! local logical for debug prints 
     208      !!---------------------------------------------------------------------- 
     209      IF( wrk_in_use(2, 1, 2, 3, 4, 5) ) THEN 
     210         CALL ctl_stop('dom_vvl_sf_nxt: requested workspace arrays unavailable')   ;   RETURN 
     211      END IF 
     212 
     213      IF(kt == nit000)   THEN 
     214         IF(lwp) WRITE(numout,*) 
     215         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' 
     216         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 
     217      ENDIF 
     218 
     219      ! ******************************* ! 
     220      ! After acale factors at t-points ! 
     221      ! ******************************* ! 
     222      !                                             !----------------------------! 
     223      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! ZTILDE or LAYER coordinate ! 
     224         !                                          !----------------------------! 
     225 
     226         ! I - Initialisations 
     227         ! =================== 
     228         IF( kt == nit000 ) THEN 
     229            ! - ML - In the future, this should be tunable by the user 
     230            IF( ln_vvl_ztilde ) THEN       ! ZTILDE CASE 
     231               ! DO jj = 1, jpj 
     232               !    DO ji = 1, jpi 
     233               !       frq_restore_hdv(ji,jj) = 2.e0 * rpi / 86400.e0 / 5.e0   & 
     234               !          &                     * MAX( SIN( ABS( gphit(ji,jj) ) / .5e0, 1.e0 / 6.e0) ) 
     235               !    END DO 
     236               ! END DO 
     237               ! frq_restore_e3t(:,:) = ( 1.e0 / 6.e0 ) * frq_restore_hdv(:,:) 
     238               frq_restore_e3t(:,:) = 2.e0 * rpi / ( 30.e0 * 86400.e0 ) 
     239               frq_restore_hdv(:,:) = 2.e0 * rpi / (  5.e0 * 86400.e0 ) 
     240               ! frq_restore_hdv(:,:) = 2.e0 * rpi / (  2.e0 * 86400.e0 ) 
     241            ELSE                           ! LAYER CASE 
     242               frq_restore_e3t(:,:) = 0.e0 
     243               frq_restore_hdv(:,:) = 0.e0 
     244            ENDIF 
     245 
     246         ENDIF 
     247 
     248         ! II - Low frequency horizontal divergence 
     249         ! ======================================== 
     250 
     251         ! 1 - barotropic divergence 
     252         ! ------------------------- 
     253         zhdiv(:,:) = 0. 
     254         zht(:,:)   = 0. 
     255         DO jk = 1, jpkm1 
     256            zhdiv(:,:) = zhdiv(:,:) + hdivn(:,:,jk) * fse3t_n(:,:,jk) 
     257            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
     258         END DO 
     259         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) ) 
     260 
     261         ! 2 - restoring equation 
     262         ! ---------------------- 
     263         IF( kt .GT. nit000 ) THEN 
     264            DO jk = 1, jpkm1 
     265               hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_restore_hdv(:,:)   & 
     266                  &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
     267            END DO 
     268         ENDIF 
     269 
     270         ! III - after z_tilde increments of vertical scale factors 
     271         ! ========================================================= 
     272         e3t_t_a(:,:,:) = 0.e0   ! e3t_t_a used to store tendency terms 
     273 
     274         ! 1 - High frequency divergence term 
     275         ! ---------------------------------- 
     276         DO jk = 1, jpkm1 
     277            e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
     278         END DO 
     279 
     280         ! 2 - Restoring term 
     281         ! ------------------ 
     282         DO jk = 1, jpk 
     283            e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - frq_restore_e3t(:,:) * e3t_t_b(:,:,jk) 
     284         END DO 
     285 
     286         ! 3 - Thickness diffusion term 
     287         ! ---------------------------- 
     288         zwu(:,:) = 0.e0 
     289         zwv(:,:) = 0.e0 
     290 
     291         ! a - first derivative: diffusive fluxes 
     292         DO jk = 1, jpkm1 
     293            DO jj = 1, jpjm1 
     294               DO ji = 1, fs_jpim1   ! vector opt. 
     295                  un_td(ji,jj,jk) = ahe3 * umask(ji,jj,jk) * e1ur(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji+1,jj  ,jk) ) 
     296                  vn_td(ji,jj,jk) = ahe3 * vmask(ji,jj,jk) * e2vr(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji  ,jj+1,jk) ) 
     297                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
     298                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
     299               END DO 
     300            END DO 
     301         END DO 
     302         ! b - correction for last oceanic u-v points 
     303         DO jj = 1, jpj 
     304            DO ji = 1, jpi 
     305               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
     306               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
     307            END DO 
     308         END DO 
     309         ! c - second derivative: divergence of diffusive fluxes 
     310         DO jk = 1, jpkm1 
     311            DO jj = 2, jpjm1 
     312               DO ji = fs_2, fs_jpim1   ! vector opt. 
     313                  e3t_t_a(ji,jj,jk) = e3t_t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
     314                     &                                      + vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk) )  & 
     315                     &                                    * e12t_1(ji,jj) 
     316               END DO 
     317            END DO 
     318         END DO 
     319         ! d - thickness diffusion equivalent transport: boundary conditions 
     320         !     (stored for tracer advaction and continuity equation) 
     321         CALL lbc_lnk( un_td , 'U' , -1.) 
     322         CALL lbc_lnk( vn_td , 'V' , -1.) 
     323 
     324         ! 4 - Time stepping of baroclinic scale factors 
     325         ! --------------------------------------------- 
     326         ! Leapfrog time stepping 
     327         ! ~~~~~~~~~~~~~~~~~~~~~~ 
     328         IF( neuler == 0 .AND. kt == nit000 ) THEN 
     329            z2dt =  rdt 
     330         ELSE 
     331            z2dt = 2.e0 * rdt 
     332         ENDIF 
     333         CALL lbc_lnk( e3t_t_a(:,:,:), 'T', 1. )    ! - ML - Do not use lbc_lnk_e3: e3t_t_a is a tendency term at this stage 
     334         e3t_t_a(:,:,:) = e3t_t_b(:,:,:) + z2dt * tmask(:,:,:) * e3t_t_a(:,:,:) 
     335         fse3t_a(:,:,:) = fse3t_0(:,:,:) + e3t_t_a(:,:,:) 
     336 
     337         ! Maximum deformation control 
     338         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
     339         ! - ML - Should perhaps be put in the namelist 
     340         z_def_max = 0.9e0 
     341         ze3t(:,:,jpk) = 0.e0 
     342         DO jk = 1, jpkm1 
     343            ze3t(:,:,jk) = e3t_t_a(:,:,jk) / fse3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 
     344         END DO 
     345         z_tmax = MAXVAL( ze3t(:,:,:) ) 
     346         z_tmin = MINVAL( ze3t(:,:,:) ) 
     347         ! - ML - test: for the moment, stop simulation for too large e3_t variations 
     348         IF( ( z_tmax .GT. z_def_max ) .OR. ( z_tmin .LT. - z_def_max ) ) THEN 
     349            ijk_max = MAXLOC( ze3t(:,:,:) ) 
     350            ijk_min = MINLOC( ze3t(:,:,:) ) 
     351            WRITE(numout, *) 'MAX( e3t_t_a(:,:,:) / fse3t_0(:,:,:) ) =', z_tmax 
     352            WRITE(numout, *) 'at i, j, k=', ijk_max 
     353            WRITE(numout, *) 'MIN( e3t_t_a(:,:,:) / fse3t_0(:,:,:) ) =', z_tmin 
     354            WRITE(numout, *) 'at i, j, k=', ijk_min             
     355            CALL ctl_stop('MAX( ABS( e3t_t_a(:,:,:) ) / fse3t_0(:,:,:) ) too high') 
     356         ENDIF 
     357         ! - ML - end test 
     358         ! - ML - This will cause a baroclinicity error if the ctl_stop above is not used 
     359         e3t_t_a(:,:,:) = MIN( e3t_t_a(:,:,:), ( 1.e0 + z_def_max ) * fse3t_0(:,:,:) ) 
     360         e3t_t_a(:,:,:) = MAX( e3t_t_a(:,:,:), ( 1.e0 - z_def_max ) * fse3t_0(:,:,:) ) 
     361 
     362         ! Boundary conditions 
     363         ! ~~~~~~~~~~~~~~~~~~~ 
     364         ! - ML - think it's not necessary 
     365         ! CALL lbc_lnk( e3t_t_a(:,:,:), 'T', 1. )    ! - ML - Do not use lbc_lnk_e3: e3t_t_a is a level thickness ANOMALY 
     366 
     367         ! IV - Barotropic repartition of the sea surface height over the baroclinic profile 
     368         ! ================================================================================= 
     369         ! add e3t(n-1) "star" Asselin-filtered 
     370         DO jk = 1, jpkm1 
     371            ! - ML - : multiplication by tmask not necessary a priori. Just to be sure for the moment. 
     372            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + ( fse3t_b(:,:,jk) - fse3t_0(:,:,jk) - e3t_t_b(:,:,jk) ) * tmask(:,:,jk) 
     373         END DO 
     374         ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n) 
     375         ! - ML - baroclinicity error should be better treated in the future 
     376         !        i.e. locally and not spread over the water column. 
     377         zht(:,:) = 0. 
     378         DO jk = 1, jpkm1 
     379            zht(:,:) = zht(:,:) + e3t_t_a(:,:,jk) * tmask(:,:,jk) 
     380         END DO 
     381         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
     382         DO jk = 1, jpkm1 
     383            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     384         END DO 
     385 
     386         !                                !------------------! 
     387      ELSEIF( ln_vvl_zstar ) THEN         ! Zstar coordinate ! 
     388         !                                !------------------! 
     389         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
     390         DO jk = 1, jpkm1 
     391            fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) 
     392         END DO 
     393 
     394      ENDIF 
     395 
     396      IF( ln_debug ) THEN   ! - ML - test: control prints for debuging 
     397         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     398            WRITE(numout, *) 'kt  =', kt 
     399            WRITE(numout, *) 'MAXVAL(abs(ht_0-SUM(e3t_t_a))) =',   & 
     400               &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) - zht(:,:) ) ) 
     401         END IF 
     402         zht(:,:) = 0.e0 
     403         DO jk = 1, jpkm1 
     404            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
     405         END DO 
     406         WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =',   & 
     407            &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 
     408         zht(:,:) = 0.e0 
     409         DO jk = 1, jpkm1 
     410            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk) 
     411         END DO 
     412         WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_a))) =',   & 
     413            &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 
     414      END IF 
     415 
     416      ! *********************************** ! 
     417      ! After scale factors at u- v- points ! 
     418      ! *********************************** ! 
     419 
     420      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' ) 
     421      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' ) 
     422 
     423      IF( wrk_not_released(2, 1, 2, 3, 4, 5) ) THEN 
     424         CALL ctl_stop( 'dom_vvl_sf_nxt: failed to release workspace arrays' ) 
     425      ENDIF 
     426 
     427   END SUBROUTINE dom_vvl_sf_nxt 
     428 
     429 
     430   SUBROUTINE dom_vvl_sf_swp( kt ) 
     431      !!---------------------------------------------------------------------- 
     432      !!                ***  ROUTINE dom_vvl_sf_swp  *** 
     433      !!                    
     434      !! ** Purpose :  compute time filter and swap of scale factors  
     435      !!               compute all depths and related variables for next time step 
     436      !!               write outputs and restart file 
     437      !! 
     438      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation 
     439      !!               - reconstruct scale factor at other grid points (interpolate) 
     440      !!               - recompute depths and water height fields 
     441      !! 
     442      !! ** Action  :  - fse3t_(b/n), e3t_t_(b/n) and fse3(u/v)_n ready for next time step 
     443      !!               - Recompute: 
     444      !!                    fse3(u/v)_b        
     445      !!                    fse3w_n            
     446      !!                    fse3(u/v)w_b       
     447      !!                    fse3(u/v)w_n       
     448      !!                    fsdept_n, fsdepw_n  and fsde3w_n 
     449      !!                    h(u/v) and h(u/v)r 
     450      !! 
     451      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     452      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling. 
     453      !!---------------------------------------------------------------------- 
     454      USE oce, ONLY:   z_e3t_def => ta                ! square of baroclinic scale factors deformation (in %) 
     455      !! * Arguments 
     456      INTEGER, INTENT( in )               :: kt       ! time step 
     457      !! * Local declarations 
     458      INTEGER                             :: jk       ! dummy loop indices 
     459      !!---------------------------------------------------------------------- 
     460 
     461      IF( kt == nit000 )   THEN 
     462         IF(lwp) WRITE(numout,*) 
     463         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors' 
     464         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step' 
     465      ENDIF 
     466 
     467      ! Time filter and swap of scale factors 
     468      ! ===================================== 
     469      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt. 
     470      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     471         IF( neuler == 0 .AND. kt == nit000 ) THEN 
     472            e3t_t_b(:,:,:) = e3t_t_n(:,:,:) 
     473         ELSE 
     474            e3t_t_b(:,:,:) = e3t_t_n(:,:,:) + atfp * (e3t_t_b(:,:,:) - 2.e0 * e3t_t_n(:,:,:) + e3t_t_a(:,:,:) ) 
     475         ENDIF 
     476         e3t_t_n(:,:,:) = e3t_t_a(:,:,:) 
     477      ENDIF 
     478      fse3t_n(:,:,:) = fse3t_a(:,:,:) 
     479      fse3u_n(:,:,:) = fse3u_a(:,:,:) 
     480      fse3v_n(:,:,:) = fse3v_a(:,:,:) 
     481 
     482      ! Compute all missing vertical scale factor and depths 
     483      ! ==================================================== 
     484      ! Horizontal scale factor interpolations 
     485      ! -------------------------------------- 
     486      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt 
     487      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  ) 
     488      ! Vertical scale factor interpolations 
     489      ! ------------------------------------ 
     490      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  ) 
     491      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 
     492      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 
     493      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 
     494      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 
     495      ! t- and w- points depth 
     496      ! ---------------------- 
     497      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1) 
     498      fsdepw_n(:,:,1) = 0.e0 
     499      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 
     500      DO jk = 2, jpk 
     501         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 
     502         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 
     503         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:) 
    167504      END DO 
    168       CALL lbc_lnk( sshu_n, 'U', 1. )   ;   CALL lbc_lnk( sshu_b, 'U', 1. )      ! lateral boundary conditions 
    169       CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. ) 
    170       CALL lbc_lnk( sshf_n, 'F', 1. ) 
    171  
    172                                                 ! initialise before scale factors at (u/v)-points 
    173       ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    174       DO jk = 1, jpkm1 
    175          DO jj = 1, jpjm1 
    176             DO ji = 1, jpim1 
    177                zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    178                zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    179                zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    180                fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    181                fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
     505      ! Local depth and Inverse of the local depth of the water column at u- and v- points 
     506      ! ---------------------------------------------------------------------------------- 
     507      hu(:,:) = 0. 
     508      hv(:,:) = 0. 
     509      DO jk = 1, jpk 
     510         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 
     511         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
     512      END DO 
     513      ! Inverse of the local depth 
     514      hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) ) 
     515      hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - umask(:,:,1) ) 
     516 
     517      ! Write outputs 
     518      ! ============= 
     519      ! - ML - add output variables in xml file for all configurations 
     520      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - fse3t_0(:,:,:) ) / fse3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
     521      CALL iom_put( "e3t_n"  , fse3t_n  (:,:,:) ) 
     522      CALL iom_put( "dept"   , fsdept_n (:,:,:) ) 
     523      CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) ) 
     524 
     525      ! write restart file 
     526      ! ================== 
     527      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 
     528 
     529   END SUBROUTINE dom_vvl_sf_swp 
     530 
     531 
     532   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 
     533      !!--------------------------------------------------------------------- 
     534      !!                  ***  ROUTINE dom_vvl__interpol  *** 
     535      !! 
     536      !! ** Purpose :   interpolate scale factors from one grid point to another 
     537      !! 
     538      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0) 
     539      !!                - horizontal interpolation: grid cell surface averaging 
     540      !!                - vertical interpolation: simple averaging 
     541      !!---------------------------------------------------------------------- 
     542      !! * Arguments 
     543      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated 
     544      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3 
     545      CHARACTER(len=1), INTENT( in )                    ::  pout       ! grid point of out scale factors 
     546      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
     547      !! * Local declarations 
     548      INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
     549      !!---------------------------------------------------------------------- 
     550      SELECT CASE ( pout ) 
     551         !               ! ------------------------------------- ! 
     552      CASE( 'U' )        ! interpolation from T-point to U-point ! 
     553         !               ! ------------------------------------- ! 
     554         ! horizontal surface weighted interpolation 
     555         DO jk = 1, jpkm1 
     556            DO jj = 2, jpjm1 
     557               DO ji = 1, fs_jpim1   ! vector opt. 
     558                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * e12u_1(ji,jj)                                          & 
     559                     &                  * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - fse3t_0(ji  ,jj,jk) )     & 
     560                     &                      + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - fse3t_0(ji+1,jj,jk) ) ) 
     561               END DO 
    182562            END DO 
    183563         END DO 
    184       END DO 
    185       CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions 
    186       CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
    187       ! Add initial scale factor to scale factor anomaly 
    188       fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
    189       fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 
     564         ! boundary conditions 
     565         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. ) 
     566         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3u_0(:,:,:) 
     567         !               ! ------------------------------------- ! 
     568      CASE( 'V' )        ! interpolation from T-point to V-point ! 
     569         !               ! ------------------------------------- ! 
     570         ! horizontal surface weighted interpolation 
     571         DO jk = 1, jpkm1 
     572            DO jj = 1, jpjm1 
     573               DO ji = fs_2, fs_jpim1   ! vector opt. 
     574                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * e12v_1(ji,jj)                                          & 
     575                     &                  * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - fse3t_0(ji,jj  ,jk) )     & 
     576                     &                      + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - fse3t_0(ji,jj+1,jk) ) ) 
     577               END DO 
     578            END DO 
     579         END DO 
     580         ! boundary conditions 
     581         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. ) 
     582         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3v_0(:,:,:) 
     583         !               ! ------------------------------------- ! 
     584      CASE( 'F' )        ! interpolation from U-point to F-point ! 
     585         !               ! ------------------------------------- ! 
     586         ! horizontal surface weighted interpolation 
     587         DO jk = 1, jpkm1 
     588            DO jj = 1, jpjm1 
     589               DO ji = 1, fs_jpim1   ! vector opt. 
     590                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * umask(ji,jj+1,jk) * e12f_1(ji,jj)                      & 
     591                     &                  * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - fse3u_0(ji,jj  ,jk) )     & 
     592                     &                      + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - fse3u_0(ji,jj+1,jk) ) ) 
     593               END DO 
     594            END DO 
     595         END DO 
     596         ! boundary conditions 
     597         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. ) 
     598         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3f_0(:,:,:) 
     599         !               ! ------------------------------------- ! 
     600      CASE( 'W' )        ! interpolation from T-point to W-point ! 
     601         !               ! ------------------------------------- ! 
     602         ! vertical simple interpolation 
     603         pe3_out(:,:,1) = fse3w_0(:,:,1) + pe3_in(:,:,1) - fse3t_0(:,:,1) 
     604         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing 
     605         DO jk = 2, jpk 
     606            pe3_out(:,:,jk) = fse3w_0(:,:,jk) + ( 1. - 0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3t_0(:,:,jk-1) )   & 
     607               &                              + (      0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3t_0(:,:,jk  ) ) 
     608         END DO 
     609         !               ! -------------------------------------- ! 
     610      CASE( 'UW' )       ! interpolation from U-point to UW-point ! 
     611         !               ! -------------------------------------- ! 
     612         ! vertical simple interpolation 
     613         pe3_out(:,:,1) = fse3uw_0(:,:,1) + pe3_in(:,:,1) - fse3u_0(:,:,1) 
     614         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing 
     615         DO jk = 2, jpk 
     616            pe3_out(:,:,jk) = fse3uw_0(:,:,jk) + ( 1. - 0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3u_0(:,:,jk-1) )   & 
     617               &                               + (      0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3u_0(:,:,jk  ) ) 
     618         END DO 
     619         !               ! -------------------------------------- ! 
     620      CASE( 'VW' )       ! interpolation from V-point to VW-point ! 
     621         !               ! -------------------------------------- ! 
     622         ! vertical simple interpolation 
     623         pe3_out(:,:,1) = fse3vw_0(:,:,1) + pe3_in(:,:,1) - fse3v_0(:,:,1) 
     624         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing 
     625         DO jk = 2, jpk 
     626            pe3_out(:,:,jk) = fse3vw_0(:,:,jk) + ( 1. - 0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3v_0(:,:,jk-1) )   & 
     627               &                               + (      0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3v_0(:,:,jk  ) ) 
     628         END DO 
     629      END SELECT 
     630 
     631   END SUBROUTINE dom_vvl_interpol 
     632 
     633 
     634   SUBROUTINE dom_vvl_rst( kt, cdrw ) 
     635      !!--------------------------------------------------------------------- 
     636      !!                   ***  ROUTINE dom_vvl_rst  *** 
     637      !!                      
     638      !! ** Purpose :   Read or write VVL file in restart file 
     639      !! 
     640      !! ** Method  :   use of IOM library 
     641      !!                if the restart does not contain vertical scale factors, 
     642      !!                they are set to the _0 values 
     643      !!                if the restart does not contain vertical scale factors increments (z_tilde), 
     644      !!                they are set to 0. 
     645      !!---------------------------------------------------------------------- 
     646      !! * Arguments 
     647      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     648      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     649      !! * Local declarations 
     650      INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
     651      !!---------------------------------------------------------------------- 
    190652      ! 
    191       IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays') 
    192       ! 
    193    END SUBROUTINE dom_vvl 
     653      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     654         !                                   ! =============== 
     655         IF( ln_rstart ) THEN                   !* Read the restart file 
     656            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. ) 
     657            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. ) 
     658            id3 = iom_varid( numror, 'e3t_t_b', ldstop = .FALSE. ) 
     659            id4 = iom_varid( numror, 'e3t_t_n', ldstop = .FALSE. ) 
     660            id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. ) 
     661 
     662            !                         ! ----------- ! 
     663            IF( ln_vvl_zstar ) THEN   ! z_star case ! 
     664               !                      ! ----------- ! 
     665               IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
     666                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
     667                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 
     668                  IF( neuler == 0 ) THEN 
     669                     fse3t_b(:,:,:) = fse3t_n(:,:,:) 
     670                  ENDIF 
     671               ELSE                                 ! one at least array is missing 
     672                  CALL ctl_stop( 'vvl cannot restart from a non vvl run' ) 
     673               ENDIF 
     674               IF( MIN( id3, id4, id5 ) > 0 ) CALL ctl_stop( 'z_star coordinate cannot restart from a z_tilde run' ) 
     675               !                      ! ------------ ! 
     676            ELSE                      ! z_tilde case ! 
     677               !                      ! ------------ ! 
     678               IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
     679                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
     680                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 
     681                  IF( neuler == 0 ) THEN 
     682                     fse3t_b(:,:,:) = fse3t_n(:,:,:) 
     683                  ENDIF 
     684               ELSE                                 ! one at least array is missing 
     685                  CALL ctl_stop( 'vvl cannot restart from a non vvl run' ) 
     686               ENDIF 
     687               IF( MIN( id3, id4, id5 ) > 0 ) THEN  ! all required arrays exist 
     688                  CALL iom_get( numror, jpdom_autoglo, 'e3t_t_b', e3t_t_b(:,:,:) ) 
     689                  CALL iom_get( numror, jpdom_autoglo, 'e3t_t_n', e3t_t_n(:,:,:) ) 
     690                  CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) ) 
     691               ELSE                                 ! one at least array is missing 
     692                  e3t_t_b(:,:,:) = 0.e0 
     693                  e3t_t_n(:,:,:) = 0.e0 
     694                  hdiv_lf(:,:,:) = 0.e0 
     695               ENDIF 
     696 
     697            ENDIF 
     698 
     699         ELSE                                   !* Initialize at "rest" 
     700            fse3t_b(:,:,:) = fse3t_0(:,:,:) 
     701            fse3t_n(:,:,:) = fse3t_0(:,:,:) 
     702            e3t_t_b(:,:,:) = 0.e0 
     703            e3t_t_n(:,:,:) = 0.e0 
     704            hdiv_lf(:,:,:) = 0.e0 
     705         ENDIF 
     706 
     707      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     708         !                                   ! =================== 
     709         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 
     710         !                                         ! ---------------------- ! 
     711         !                                         ! z_star & z_tilde cases ! 
     712         !                                         ! ---------------------- ! 
     713         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 
     714         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) ) 
     715         !                                           ! ----------------- ! 
     716         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde case only ! 
     717            !                                        ! ----------------- ! 
     718            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_b', e3t_t_b(:,:,:) ) 
     719            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_n', e3t_t_n(:,:,:) ) 
     720            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 
     721         ENDIF 
     722 
     723      ENDIF 
     724 
     725   END SUBROUTINE dom_vvl_rst 
     726 
     727 
     728   SUBROUTINE dom_vvl_ctl 
     729      !!--------------------------------------------------------------------- 
     730      !!                  ***  ROUTINE dom_vvl_ctl  *** 
     731      !!                 
     732      !! ** Purpose :   Control the consistency between namelist options for  
     733      !!                vertical coordinate and set nvvl 
     734      !!---------------------------------------------------------------------- 
     735      INTEGER ::   ioptio 
     736 
     737      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ahe3! , ln_vvl_kepe 
     738      !!----------------------------------------------------------------------  
     739 
     740      REWIND ( numnam )               ! Read Namelist nam_vvl : vertical coordinate 
     741      READ   ( numnam, nam_vvl ) 
     742 
     743      IF(lwp) THEN                    ! Namelist print 
     744         WRITE(numout,*) 
     745         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate' 
     746         WRITE(numout,*) '~~~~~~~~~~~' 
     747         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate' 
     748         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar 
     749         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde 
     750         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer 
     751         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation' 
     752         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe 
     753         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient' 
     754         WRITE(numout,*) '                                         ahe3           = ', ahe3 
     755      ENDIF 
     756 
     757      ioptio = 0                      ! Parameter control 
     758      IF( ln_vvl_zstar     )   ioptio = ioptio + 1 
     759      IF( ln_vvl_ztilde    )   ioptio = ioptio + 1 
     760      IF( ln_vvl_layer     )   ioptio = ioptio + 1 
     761 
     762      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
     763 
     764      IF( ln_vvl_zstar  )   nvvl = 1 
     765      IF( ln_vvl_ztilde )   nvvl = 2 
     766      IF( ln_vvl_layer  )   nvvl = 3 
     767 
     768      IF(lwp) THEN                   ! Print the choice 
     769         WRITE(numout,*) 
     770         IF( nvvl ==  1        ) WRITE(numout,*) '              zstar vertical coordinate is used' 
     771         IF( nvvl ==  2        ) WRITE(numout,*) '              ztilde vertical coordinate is used' 
     772         IF( nvvl ==  3        ) WRITE(numout,*) '              layer vertical coordinate is used' 
     773         ! - ML - Option not developed yet 
     774         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used' 
     775         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used' 
     776      ENDIF 
     777 
     778   END SUBROUTINE dom_vvl_ctl 
     779 
    194780 
    195781#else 
     782 
     783 
    196784   !!---------------------------------------------------------------------- 
    197785   !!   Default option :                                      Empty routine 
    198786   !!---------------------------------------------------------------------- 
    199 CONTAINS 
    200    SUBROUTINE dom_vvl 
    201    END SUBROUTINE dom_vvl 
     787   SUBROUTINE dom_vvl_init 
     788      WRITE(*,*) 'dom_vvl_init: You should not have seen this print! error?' 
     789   END SUBROUTINE dom_vvl_init 
     790   SUBROUTINE dom_vvl_sf_nxt( kt )  
     791      !! * Arguments 
     792      INTEGER, INTENT( in )  ::    kt 
     793      WRITE(*,*) 'dom_vvl_sf_nxt: You should not have seen this print! error?', kt 
     794   END SUBROUTINE dom_vvl_sf_nxt 
     795   SUBROUTINE dom_vvl_sf_swp( kt )  
     796      !! * Arguments 
     797      INTEGER, INTENT( in )  ::    kt 
     798      WRITE(*,*) 'dom_vvl_sf_swp: You should not have seen this print! error?', kt 
     799   END SUBROUTINE dom_vvl_sf_swp 
     800   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 
     801      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in 
     802      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out 
     803      CHARACTER(len=1), INTENT( in )                    ::  pout 
     804      WRITE(*,*) 'dom_vvl_interpol: You should not have seen this print! error?' 
     805   END SUBROUTINE dom_vvl_interpol 
     806 
     807 
    202808#endif 
    203809 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90

    r2528 r2905  
    1919#   define  fse3uw_0(i,j,k)  e3uw(i,j,k) 
    2020#   define  fse3vw_0(i,j,k)  e3vw(i,j,k) 
     21 
    2122#if defined key_vvl 
    22 ! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._1) 
    23 #   define  fsdept(i,j,k)  gdept_1(i,j,k) 
    24 #   define  fsdepw(i,j,k)  gdepw_1(i,j,k) 
    25 #   define  fsde3w(i,j,k)  gdep3w_1(i,j,k) 
    26 #   define  fse3t(i,j,k)   e3t_1(i,j,k) 
    27 #   define  fse3u(i,j,k)   e3u_1(i,j,k) 
    28 #   define  fse3v(i,j,k)   e3v_1(i,j,k) 
    29 #   define  fse3f(i,j,k)   e3f_1(i,j,k) 
    30 #   define  fse3w(i,j,k)   e3w_1(i,j,k) 
    31 #   define  fse3uw(i,j,k)  e3uw_1(i,j,k) 
    32 #   define  fse3vw(i,j,k)  e3vw_1(i,j,k) 
     23! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._n) 
    3324 
    3425#   define  fse3t_b(i,j,k)   e3t_b(i,j,k) 
    3526#   define  fse3u_b(i,j,k)   e3u_b(i,j,k) 
    3627#   define  fse3v_b(i,j,k)   e3v_b(i,j,k) 
    37 #   define  fse3uw_b(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k))) 
    38 #   define  fse3vw_b(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k))) 
     28#   define  fse3uw_b(i,j,k)  e3uw_b(i,j,k) 
     29#   define  fse3vw_b(i,j,k)  e3vw_b(i,j,k) 
    3930 
    40 #   define  fsdept_n(i,j,k)  (fsdept_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k))) 
    41 #   define  fsdepw_n(i,j,k)  (fsdepw_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k))) 
    42 #   define  fsde3w_n(i,j,k)  (fsde3w_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k))-sshn(i,j)) 
    43 #   define  fse3t_n(i,j,k)   (fse3t_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k))) 
    44 #   define  fse3u_n(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_n(i,j)*muu(i,j,k))) 
    45 #   define  fse3v_n(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k))) 
    46 #   define  fse3f_n(i,j,k)   (fse3f_0(i,j,k)*(1.+sshf_n(i,j)*muf(i,j,k))) 
    47 #   define  fse3w_n(i,j,k)   (fse3w_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k))) 
    48 #   define  fse3uw_n(i,j,k)  (fse3uw_0(i,j,k)*(1.+sshu_n(i,j)*muu(i,j,k))) 
    49 #   define  fse3vw_n(i,j,k)  (fse3vw_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k))) 
     31#   define  fsdept_n(i,j,k)  gdept_n(i,j,k) 
     32#   define  fsdepw_n(i,j,k)  gdepw_n(i,j,k) 
     33#   define  fsde3w_n(i,j,k)  gdep3w_n(i,j,k) 
     34#   define  fse3t_n(i,j,k)   e3t_n(i,j,k) 
     35#   define  fse3u_n(i,j,k)   e3u_n(i,j,k) 
     36#   define  fse3v_n(i,j,k)   e3v_n(i,j,k) 
     37#   define  fse3f_n(i,j,k)   e3f_n(i,j,k) 
     38#   define  fse3w_n(i,j,k)   e3w_n(i,j,k) 
     39#   define  fse3uw_n(i,j,k)  e3uw_n(i,j,k) 
     40#   define  fse3vw_n(i,j,k)  e3vw_n(i,j,k) 
    5041 
    51 #   define  fse3t_m(i,j,k)   (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k))) 
     42#   define  fse3t_a(i,j,k)   e3t_a(i,j,k) 
     43#   define  fse3u_a(i,j,k)   e3u_a(i,j,k) 
     44#   define  fse3v_a(i,j,k)   e3v_a(i,j,k) 
    5245 
    53 #   define  fse3t_a(i,j,k)   (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 
    54 #   define  fse3u_a(i,j,k)   (fse3u_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k))) 
    55 #   define  fse3v_a(i,j,k)   (fse3v_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k))) 
     46#   define  fse3t_m(i,j)   e3t_m(i,j) 
     47 
     48! This part should be removed one day ... 
     49#   define  fsdept(i,j,k)  fsdept_n(i,j,k) 
     50#   define  fsdepw(i,j,k)  fsdepw_n(i,j,k) 
     51#   define  fsde3w(i,j,k)  fsde3w_n(i,j,k) 
     52#   define  fse3t(i,j,k)   fse3t_n(i,j,k)  
     53#   define  fse3u(i,j,k)   fse3u_n(i,j,k)  
     54#   define  fse3v(i,j,k)   fse3v_n(i,j,k)  
     55#   define  fse3f(i,j,k)   fse3f_n(i,j,k)  
     56#   define  fse3w(i,j,k)   fse3w_n(i,j,k)  
     57#   define  fse3uw(i,j,k)  fse3uw_n(i,j,k) 
     58#   define  fse3vw(i,j,k)  fse3vw_n(i,j,k) 
    5659 
    5760#else 
    … …  
    8588#   define  fse3vw_n(i,j,k)  fse3vw_0(i,j,k) 
    8689 
    87 #   define  fse3t_m(i,j,k)   fse3t_0(i,j,k) 
     90#   define  fse3t_m(i,j)     fse3t_0(i,j,1) 
    8891 
    8992#   define  fse3t_a(i,j,k)   fse3t_0(i,j,k) 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r2723 r2905  
    9292      !!               un,vn   now horizontal velocity of next time-step 
    9393      !!---------------------------------------------------------------------- 
    94       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    9594      USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace 
    96       USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3 
    9795      ! 
    9896      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    … …  
    104102      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars 
    105103      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      - 
    106       REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1 
     104      REAL(wp) ::   zec                         !   -      - 
    107105      !!---------------------------------------------------------------------- 
    108  
    109       IF( wrk_in_use(2, 1,2,3) ) THEN 
    110          CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN 
    111       ENDIF 
    112106 
    113107      IF( kt == nit000 ) THEN 
    … …  
    185179         hvr_e(:,:) = hvr(:,:) 
    186180         DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    187             ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    188             va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     181            ua_e(:,:) = ua_e(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     182            va_e(:,:) = va_e(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    189183         END DO 
    190184         ua_e(:,:) = ua_e(:,:) * hur(:,:) 
    … …  
    239233            !                             ! ================! 
    240234            ! Before scale factor at t-points 
    241             ! ------------------------------- 
     235            ! (used as a now filtered scale factor until the swap) 
     236            ! ---------------------------------------------------- 
    242237            DO jk = 1, jpkm1 
    243238               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   & 
    … …  
    245240                  &                         - 2.e0 * fse3t_n(:,:,jk)            ) 
    246241            ENDDO 
    247             ! Add volume filter correction only at the first level of t-point scale factors 
     242            ! Add volume filter correction: comatibility with tracer advection scheme 
     243            ! => time filter + conservation correction (only at the first level) 
    248244            zec = atfp * rdt / rau0 
    249245            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    250             ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 
    251             zs_t  (:,:) =       e1t(:,:) * e2t(:,:) 
    252             zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:) 
    253             zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:) 
     246            ! swap of emp fields 
     247            emp_b(:,:) = emp(:,:) 
    254248            ! 
    255249            IF( ln_dynadv_vec ) THEN 
    256250               ! Before scale factor at (u/v)-points 
    257251               ! ----------------------------------- 
    258                ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    259                DO jk = 1, jpkm1 
    260                   DO jj = 1, jpjm1 
    261                      DO ji = 1, jpim1 
    262                         zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    263                         zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    264                         zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    265                         fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    266                         fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
    267                      END DO 
    268                   END DO 
    269                END DO 
    270                ! lateral boundary conditions 
    271                CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) 
    272                CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
    273                ! Add initial scale factor to scale factor anomaly 
    274                fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
    275                fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 
     252               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 
     253               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 
    276254               ! Leap-Frog - Asselin filter and swap: applied on velocity 
    277255               ! ----------------------------------- 
    … …  
    291269               ! 
    292270            ELSE 
    293                ! Temporary filered scale factor at (u/v)-points (will become before scale factor) 
    294                !----------------------------------------------- 
    295                ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    296                DO jk = 1, jpkm1 
    297                   DO jj = 1, jpjm1 
    298                      DO ji = 1, jpim1 
    299                         zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    300                         zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    301                         zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    302                         ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    303                         ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
    304                      END DO 
    305                   END DO 
    306                END DO 
    307                ! lateral boundary conditions 
    308                CALL lbc_lnk( ze3u_f, 'U', 1. ) 
    309                CALL lbc_lnk( ze3v_f, 'V', 1. ) 
    310                ! Add initial scale factor to scale factor anomaly 
    311                ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 
    312                ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 
     271               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor) 
     272               !------------------------------------------------ 
     273               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' ) 
     274               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' ) 
    313275               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 
    314276               ! -----------------------------------             =========================== 
    … …  
    345307      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    346308         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    347       !  
    348       IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays') 
     309      ! 
    349310      ! 
    350311   END SUBROUTINE dyn_nxt 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r2724 r2905  
    114114      USE wrk_nemo, ONLY:   zun => wrk_2d_14 , zun_e => wrk_2d_15 , zub_e => wrk_2d_16 , zu_sum => wrk_2d_17 
    115115      USE wrk_nemo, ONLY:   zvn => wrk_2d_18 , zvn_e => wrk_2d_19 , zvb_e => wrk_2d_20 , zv_sum => wrk_2d_21 
     116      USE wrk_nemo, ONLY:   zhu_b => wrk_2d_22 , zhv_b => wrk_2d_23 
    116117      ! 
    117118      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    … …  
    127128 
    128129      IF( wrk_in_use(2,  1, 2, 3, 4, 5, 6, 7, 8, 9,10,     & 
    129          &              11,12,13,14,15,16,17,18,19,20,21 ) ) THEN 
     130         &              11,12,13,14,15,16,17,18,19,20,21,22,23 ) ) THEN 
    130131         CALL ctl_stop( 'dyn_spg_ts: requested workspace arrays unavailable' )   ;   RETURN 
    131132      ENDIF 
    … …  
    188189#endif 
    189190               !                                                                              ! now trend 
    190                zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
    191                zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
     191               zua(ji,jj) = zua(ji,jj) + fse3u_n(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
     192               zva(ji,jj) = zva(ji,jj) + fse3v_n(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    192193               !                                                                              ! now velocity  
    193                zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk) 
    194                zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)                
     194               zun(ji,jj) = zun(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) 
     195               zvn(ji,jj) = zvn(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk)                
    195196               ! 
    196 #if defined key_vvl 
    197                ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk)  
    198                vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)    
    199 #else 
    200                ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
    201                vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk) 
    202 #endif 
     197               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
     198               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk) 
    203199            END DO 
    204200         END DO 
    205201      END DO 
     202 
     203      ! before inverse water column height at u- and v- points 
     204      IF( lk_vvl ) THEN 
     205         zhu_b(:,:) = 0. 
     206         zhv_b(:,:) = 0. 
     207         DO jk = 1, jpk 
     208            zhu_b(:,:) = zhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 
     209            zhv_b(:,:) = zhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
     210         END DO 
     211         zhu_b(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1) ) 
     212         zhv_b(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1) ) 
     213      ELSE 
     214         zhu_b(:,:) = hur(:,:) 
     215         zhv_b(:,:) = hvr(:,:) 
     216      ENDIF 
    206217 
    207218      !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
    … …  
    289300            !RBbug for vvl and external mode we may need to use varying fse3 
    290301            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 
    291             zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
    292             zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
     302            zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u_n(ji,jj,mbku(ji,jj)) * zcoef  ) 
     303            zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v_n(ji,jj,mbkv(ji,jj)) * zcoef  ) 
    293304         END DO 
    294305      END DO 
    295  
    296       IF( lk_vvl ) THEN 
    297          DO jj = 2, jpjm1 
    298             DO ji = fs_2, fs_jpim1   ! vector opt. 
    299                zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
    300                   &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 
    301                zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
    302                   &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 
    303             END DO 
    304          END DO 
    305       ELSE 
    306          DO jj = 2, jpjm1 
    307             DO ji = fs_2, fs_jpim1   ! vector opt. 
    308                zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
    309                zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
    310             END DO 
    311          END DO 
    312       ENDIF 
     306       
     307      DO jj = 2, jpjm1 
     308         DO ji = fs_2, fs_jpim1   ! vector opt. 
     309            zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * zhu_b(ji,jj) 
     310            zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * zhv_b(ji,jj) 
     311         END DO 
     312      END DO 
    313313 
    314314      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
    … …  
    317317      zva(:,:) = zva(:,:) * hvr(:,:) 
    318318      ! 
    319       IF( lk_vvl ) THEN 
    320          ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 
    321          vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 
    322       ELSE 
    323          ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    324          vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
    325       ENDIF 
     319      ub_b(:,:) = ub_b(:,:) * zhu_b(:,:) 
     320      vb_b(:,:) = vb_b(:,:) * zhv_b(:,:) 
    326321 
    327322      ! ----------------------------------------------------------------------- 
    … …  
    640635            ENDIF 
    641636 
    642             IF( lk_vvl ) THEN 
    643                ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 
    644                vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 
    645             ELSE 
    646                ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    647                vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
    648             ENDIF 
    649          ELSE                                 ! neuler==0 
    650             ub_b (:,:) = un_b (:,:) 
    651             vb_b (:,:) = vn_b (:,:) 
    652          ENDIF 
     637            ub_b(:,:) = ub_b(:,:) * zhu_b(:,:) 
     638            vb_b(:,:) = vb_b(:,:) * zhv_b(:,:) 
    653639 
    654640         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r2715 r2905  
    88   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    99   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
     10   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work 
    1011   !!---------------------------------------------------------------------- 
    1112 
    … …  
    4344   PRIVATE 
    4445 
    45    PUBLIC   ssh_wzv    ! called by step.F90 
    4646   PUBLIC   ssh_nxt    ! called by step.F90 
     47   PUBLIC   wzv        ! called by step.F90 
     48   PUBLIC   ssh_swp    ! called by step.F90 
    4749 
    4850   !! * Substitutions 
    … …  
    5658CONTAINS 
    5759 
    58    SUBROUTINE ssh_wzv( kt )  
    59       !!---------------------------------------------------------------------- 
    60       !!                ***  ROUTINE ssh_wzv  *** 
     60   SUBROUTINE ssh_nxt( kt ) 
     61      !!---------------------------------------------------------------------- 
     62      !!                ***  ROUTINE ssh_nxt  *** 
    6163      !!                    
    62       !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity 
    63       !!              and update the now vertical coordinate (lk_vvl=T). 
    64       !! 
    65       !! ** Method  : - Using the incompressibility hypothesis, the vertical  
    66       !!      velocity is computed by integrating the horizontal divergence   
    67       !!      from the bottom to the surface minus the scale factor evolution. 
    68       !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     64      !! ** Purpose :   compute the after ssh (ssha) 
     65      !! 
     66      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment 
     67      !!      is computed by integrating the horizontal divergence and multiply by 
     68      !!      by the time step. 
    6969      !! 
    7070      !! ** action  :   ssha    : after sea surface height 
    71       !!                wn      : now vertical velocity 
    72       !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T) 
    73       !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points 
    7471      !! 
    7572      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    7673      !!---------------------------------------------------------------------- 
    7774      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    78       USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace 
    79       USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1 , z2d => wrk_2d_2   ! 2D workspace 
    80       ! 
    81       INTEGER, INTENT(in) ::   kt   ! time step 
    82       ! 
    83       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    84       REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars 
    85       !!---------------------------------------------------------------------- 
    86  
    87       IF( wrk_in_use(2, 1,2) ) THEN 
    88          CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN 
     75      USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1                     ! 2D workspace 
     76      ! 
     77      INTEGER, INTENT(in) ::   kt                      ! time step 
     78      !  
     79      INTEGER             ::   jk                      ! dummy loop indice 
     80      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
     81      !!---------------------------------------------------------------------- 
     82 
     83      IF( wrk_in_use(2, 1) ) THEN 
     84         CALL ctl_stop('ssh_nxt: requested workspace arrays unavailable')   ;   RETURN 
    8985      ENDIF 
    9086 
    … …  
    9288         ! 
    9389         IF(lwp) WRITE(numout,*) 
    94          IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' 
     90         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 
    9591         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    9692         ! 
    97          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
    98          ! 
    99          IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only) 
    100             DO jj = 1, jpjm1 
    101                DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible 
    102                   zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
    103                   zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
    104                   zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
    105                   sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
    106                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
    107                   sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    108                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    109                   sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    110                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
    111                   sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
    112                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
    113                END DO 
    114             END DO 
    115             CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
    116             CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
    117             DO jj = 1, jpjm1 
    118                DO ji = 1, jpim1      ! NO Vector Opt. 
    119                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
    120                        &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    121                        &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    122                        &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    123                END DO 
    124             END DO 
    125             CALL lbc_lnk( sshf_n, 'F', 1. ) 
    126          ENDIF 
    127          ! 
    128       ENDIF 
    129  
    130       !                                           !------------------------------------------! 
    131       IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case) 
    132          !                                        !------------------------------------------! 
    133          DO jk = 1, jpkm1 
    134             fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays 
    135             fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 
    136             fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 
    137             ! 
    138             fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays 
    139             fse3u (:,:,jk) = fse3u_n (:,:,jk) 
    140             fse3v (:,:,jk) = fse3v_n (:,:,jk) 
    141             fse3f (:,:,jk) = fse3f_n (:,:,jk) 
    142             fse3w (:,:,jk) = fse3w_n (:,:,jk) 
    143             fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 
    144             fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 
    145          END DO 
    146          ! 
    147          hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points) 
    148          hv(:,:) = hv_0(:,:) + sshv_n(:,:) 
    149          !                                            ! now masked inverse of the ocean depth (at u- and v-points) 
    150          hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) ) 
    151          hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) ) 
    152          !  
    15393      ENDIF 
    15494      ! 
    … …  
    184124      CALL lbc_lnk( ssha, 'T', 1. )  
    185125#endif 
    186  
    187       !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only) 
    188       IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
    189          DO jj = 1, jpjm1 
    190             DO ji = 1, jpim1      ! NO Vector Opt. 
    191                sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
    192                   &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     & 
    193                   &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
    194                sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
    195                   &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
    196                   &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
    197             END DO 
    198          END DO 
    199          CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions 
    200       ENDIF 
    201        
    202126#if defined key_asminc 
    203127      !                                                ! Include the IAU weighted SSH increment 
    … …  
    209133 
    210134      !                                           !------------------------------! 
     135      !                                           !           outputs            ! 
     136      !                                           !------------------------------! 
     137      CALL iom_put( "ssh" , sshn                  )   ! sea surface height 
     138      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
     139      ! 
     140      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
     141      ! 
     142      IF( wrk_not_released(2, 1) )   CALL ctl_stop('ssh_nxt: failed to release workspace arrays') 
     143      ! 
     144   END SUBROUTINE ssh_nxt 
     145 
     146    
     147   SUBROUTINE wzv( kt ) 
     148      !!---------------------------------------------------------------------- 
     149      !!                ***  ROUTINE wzv  *** 
     150      !!                    
     151      !! ** Purpose :   compute the now vertical velocity 
     152      !! 
     153      !! ** Method  : - Using the incompressibility hypothesis, the vertical  
     154      !!      velocity is computed by integrating the horizontal divergence   
     155      !!      from the bottom to the surface minus the scale factor evolution. 
     156      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     157      !! 
     158      !! ** action  :   wn      : now vertical velocity 
     159      !! 
     160      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     161      !!---------------------------------------------------------------------- 
     162      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
     163      USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace 
     164      USE wrk_nemo, ONLY:   zhdiv => wrk_3d_1 , z2d => wrk_2d_1   ! 2D workspace 
     165      ! 
     166      INTEGER, INTENT(in) ::   kt           ! time step 
     167      ! 
     168      INTEGER             ::   ji, jj, jk   ! dummy loop indices 
     169      REAL(wp)            ::   z1_2dt       ! local scalars 
     170      !!---------------------------------------------------------------------- 
     171       
     172      IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 1) ) THEN 
     173         CALL ctl_stop('wzv: requested workspace arrays unavailable')   ;   RETURN 
     174      ENDIF 
     175 
     176      IF( kt == nit000 ) THEN 
     177         ! 
     178         IF(lwp) WRITE(numout,*) 
     179         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 
     180         IF(lwp) WRITE(numout,*) '~~~ ' 
     181         ! 
     182         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
     183         ! 
     184      ENDIF 
     185      !                                           !------------------------------! 
    211186      !                                           !     Now Vertical Velocity    ! 
    212187      !                                           !------------------------------! 
    213       z1_2dt = 1.e0 / z2dt 
    214       DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    215          ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
    216          wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    217             &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
    218             &                         * tmask(:,:,jk) * z1_2dt 
     188      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog) 
     189      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt 
     190      ! 
     191      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
     192         DO jk = 1, jpkm1 
     193            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
     194            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
     195            DO jj = 2, jpjm1 
     196               DO ji = fs_2, fs_jpim1   ! vector opt. 
     197                  zhdiv(ji,jj,jk) = e12t_1(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
     198               END DO 
     199            END DO 
     200         END DO 
     201         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
     202         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells? 
     203         !                             ! Same question holds for hdivn. Perhaps just for security 
     204         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     205            ! computation of w 
     206            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    & 
     207               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
     208         END DO 
     209         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
     210      ELSE   ! z_star and linear free surface cases 
     211         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     212            ! computation of w 
     213            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   & 
     214               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
     215         END DO 
     216      ENDIF 
     217 
    219218#if defined key_bdy 
    220219         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
    221220#endif 
    222       END DO 
    223221 
    224222      !                                           !------------------------------! 
    225223      !                                           !           outputs            ! 
    226224      !                                           !------------------------------! 
    227       CALL iom_put( "woce", wn                    )   ! vertical velocity 
    228       CALL iom_put( "ssh" , sshn                  )   ! sea surface height 
    229       CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
     225      CALL iom_put( "woce", wn )   ! vertical velocity 
    230226      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value 
    231227         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    … …  
    237233         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    238234      ENDIF 
    239       ! 
    240       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
    241       ! 
    242       IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays') 
    243       ! 
    244    END SUBROUTINE ssh_wzv 
    245  
    246  
    247    SUBROUTINE ssh_nxt( kt ) 
     235 
     236      IF( wrk_not_released(2, 1) .OR. wrk_not_released(3, 1) )   CALL ctl_stop('wzv: failed to release workspace arrays') 
     237 
     238   END SUBROUTINE wzv 
     239 
     240 
     241   SUBROUTINE ssh_swp( kt ) 
    248242      !!---------------------------------------------------------------------- 
    249243      !!                    ***  ROUTINE ssh_nxt  *** 
    … …  
    251245      !! ** Purpose :   achieve the sea surface  height time stepping by  
    252246      !!              applying Asselin time filter and swapping the arrays 
    253       !!              ssha  already computed in ssh_wzv   
     247      !!              ssha  already computed in ssh_nxt   
    254248      !! 
    255249      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing 
    … …  
    266260      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    267261      !! 
    268       INTEGER  ::   ji, jj   ! dummy loop indices 
    269       REAL(wp) ::   zec      ! temporary scalar 
     262      REAL(wp)            ::   zec  ! temporary scalar 
    270263      !!---------------------------------------------------------------------- 
    271264 
    272265      IF( kt == nit000 ) THEN 
    273266         IF(lwp) WRITE(numout,*) 
    274          IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' 
     267         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
    275268         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    276269      ENDIF 
    277270 
    278       !                       !--------------------------! 
    279       IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points) 
    280          !                    !--------------------------! 
    281          ! 
    282          IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
    283             sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now) 
    284             sshu_n(:,:) = sshu_a(:,:) 
    285             sshv_n(:,:) = sshv_a(:,:) 
    286             DO jj = 1, jpjm1                                ! ssh now at f-point 
    287                DO ji = 1, jpim1      ! NO Vector Opt. 
    288                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 & 
    289                      &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    290                      &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    291                      &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    292                END DO 
    293             END DO 
    294             CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions 
    295             ! 
    296          ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap 
     271      IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
     272         sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now) 
     273      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap 
     274         IF( lk_vvl ) THEN                               ! before <-- now filtered 
    297275            zec = atfp * rdt / rau0 
    298             DO jj = 1, jpj 
    299                DO ji = 1, jpi                               ! before <-- now filtered 
    300                   sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   & 
    301                      &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 
    302                   sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after 
    303                   sshu_n(ji,jj) = sshu_a(ji,jj) 
    304                   sshv_n(ji,jj) = sshv_a(ji,jj) 
    305                END DO 
    306             END DO 
    307             DO jj = 1, jpjm1                                ! ssh now at f-point 
    308                DO ji = 1, jpim1      ! NO Vector Opt. 
    309                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 & 
    310                      &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    311                      &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    312                      &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    313                END DO 
    314             END DO 
    315             CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions 
    316             ! 
    317             DO jj = 1, jpjm1                                ! ssh before at u- & v-points 
    318                DO ji = 1, jpim1      ! NO Vector Opt. 
    319                   sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
    320                      &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
    321                      &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
    322                   sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
    323                      &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    324                      &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    325                END DO 
    326             END DO 
    327             CALL lbc_lnk( sshu_b, 'U', 1. ) 
    328             CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions 
    329             ! 
     276            sshb  (:,:) = sshn  (:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )   & 
     277               &                          - zec  * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     278         ELSE 
     279            sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 
    330280         ENDIF 
    331          !                    !--------------------------! 
    332       ELSE                    !        fixed levels      !     (ssh at t-point only) 
    333          !                    !--------------------------! 
    334          ! 
    335          IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
    336             sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now) 
    337             ! 
    338          ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap 
    339             DO jj = 1, jpj 
    340                DO ji = 1, jpi                               ! before <-- now filtered 
    341                   sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 
    342                   sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after 
    343                END DO 
    344             END DO 
    345          ENDIF 
    346          ! 
     281         sshn(:,:) = ssha(:,:)                           ! now <-- after 
    347282      ENDIF 
    348283      ! 
    … …  
    354289      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
    355290      ! 
    356    END SUBROUTINE ssh_nxt 
     291   END SUBROUTINE ssh_swp 
    357292 
    358293   !!====================================================================== 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r2528 r2905  
    2323   USE eosbn2          ! equation of state            (eos bn2 routine) 
    2424   USE trdmld_oce      ! ocean active mixed layer tracers trends variables 
    25    USE domvvl          ! variable volume 
    2625   USE traswp          ! swap from 4D T-S to 3D T & S and vice versa 
    2726 
    … …  
    122121                     CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb     ) 
    123122                     CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb      ) 
    124       IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 
    125123                     ! 
    126124                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , un        )     ! now fields 
    … …  
    191189                     CALL iom_get( numror, jpdom_autoglo, 'hdivb'  , hdivb   ) 
    192190                     CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    ) 
    193       IF( lk_vvl )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    194191                     ! 
    195192                     CALL iom_get( numror, jpdom_autoglo, 'un'     , un      )   ! now    fields 
    … …  
    218215         hdivb(:,:,:) = hdivn(:,:,:) 
    219216         sshb (:,:)   = sshn (:,:) 
    220          IF( lk_vvl ) THEN 
    221             DO jk = 1, jpk 
    222                fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
    223             END DO 
    224          ENDIF 
    225217      ENDIF 
    226218      ! 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r2715 r2905  
    8282   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sss_m     !: mean (nn_fsbc time-step) surface sea salinity            [psu] 
    8383   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssh_m     !: mean (nn_fsbc time-step) sea surface height                [m] 
     84   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface height                [m] 
    8485 
    8586   !! * Substitutions 
    … …  
    9697      !!                  ***  FUNCTION sbc_oce_alloc  *** 
    9798      !!--------------------------------------------------------------------- 
    98       INTEGER :: ierr(4) 
     99      INTEGER :: ierr(5) 
    99100      !!--------------------------------------------------------------------- 
    100101      ierr(:) = 0 
    … …  
    117118         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) ,                       & 
    118119         &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) ) 
     120         ! 
     121#if defined key_vvl 
     122      ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 
     123#endif 
    119124         ! 
    120125      sbc_oce_alloc = MAXVAL( ierr ) 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r2715 r2905  
    7070         ELSE                    ;   ssh_m(:,:) = sshn(:,:) 
    7171         ENDIF 
    72  
     72         ! 
     73         IF( lk_vvl )   fse3t_m(:,:) = fse3t_n(:,:,1) 
    7374         ! 
    7475      ELSE 
    … …  
    8687               CALL iom_get( numror, jpdom_autoglo, 'sss_m'  , sss_m  )   !   "         "    salinity    (T-point) 
    8788               CALL iom_get( numror, jpdom_autoglo, 'ssh_m'  , ssh_m  )   !   "         "    height      (T-point) 
     89               IF( lk_vvl ) THEN 
     90                  CALL iom_get( numror, jpdom_autoglo, 'fse3t_m', fse3t_m(:,:) ) 
     91                  !                                                       !   "         "    vertical scale factor (T-point) 
     92               ENDIF 
    8893               ! 
    8994               IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs 
    … …  
    96101                  sss_m(:,:) = zcoef * sss_m(:,:) 
    97102                  ssh_m(:,:) = zcoef * ssh_m(:,:) 
     103                  IF( lk_vvl )   fse3t_m(:,:) = zcoef * fse3t_m(:,:) 
    98104               ELSE 
    99105                  IF(lwp) WRITE(numout,*) '~~~~~~~   mean fields read in the ocean restart file' 
    … …  
    110116               ELSE                    ;   ssh_m(:,:) = zcoef *   sshn(:,:) 
    111117               ENDIF 
    112  
     118               IF( lk_vvl )   fse3t_m(:,:) = zcoef * fse3t_n(:,:,1) 
     119               ! 
    113120            ENDIF 
    114121            !                                             ! ---------------------------------------- ! 
    … …  
    120127            sss_m(:,:) = 0.e0 
    121128            ssh_m(:,:) = 0.e0 
     129            IF( lk_vvl )   fse3t_m(:,:) = 0.e0 
    122130         ENDIF 
    123131         !                                                ! ---------------------------------------- ! 
    … …  
    132140         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) 
    133141         ENDIF 
     142         IF( lk_vvl )   fse3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1) 
    134143 
    135144         !                                                ! ---------------------------------------- ! 
    … …  
    142151            ssv_m(:,:) = ssv_m(:,:) * zcoef           ! 
    143152            ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m] 
     153            IF( lk_vvl )   fse3t_m(:,:) = fse3t_m(:,:) * zcoef   ! mean vertical scale factor [m] 
    144154            ! 
    145155         ENDIF 
    … …  
    158168            CALL iom_rstput( kt, nitrst, numrow, 'sss_m'  , sss_m  ) 
    159169            CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  ) 
     170            IF( lk_vvl ) THEN 
     171               CALL iom_rstput( kt, nitrst, numrow, 'fse3t_m'  , fse3t_m(:,:)  ) 
     172            END IF 
    160173            ! 
    161174         ENDIF 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r2715 r2905  
    1414   USE oce             ! ocean dynamics and active tracers 
    1515   USE dom_oce         ! ocean space and time domain 
     16   USE domvvl          ! variable vertical scale factors 
    1617   USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine) 
    1718   USE traadv_tvd      ! TVD      scheme           (tra_adv_tvd    routine) 
    … …  
    8889         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk) 
    8990      END DO 
     91      ! 
     92      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     93         zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 
     94         zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 
     95      ENDIF 
     96      ! 
    9097      zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    9198      zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90

    r2715 r2905  
    107107            DO jj = 1, jpjm1 
    108108               DO ji = 1, fs_jpim1   ! vector opt. 
    109                   zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 
    110                   zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 
     109                  zeeu(ji,jj) = e1ur(ji,jj) * fse3u_n(ji,jj,jk) * umask(ji,jj,jk) 
     110                  zeev(ji,jj) = e2vr(ji,jj) * fse3v_n(ji,jj,jk) * vmask(ji,jj,jk) 
    111111               END DO 
    112112            END DO 
    … …  
    130130            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient 
    131131               DO ji = fs_2, fs_jpim1   ! vector opt. 
    132                   zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     132                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    133133                  zlt(ji,jj) = fsahtt(ji,jj,jk) * zbtr * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
    134134                     &                                     + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   ) 
    … …  
    148148               DO ji = fs_2, fs_jpim1   ! vector opt. 
    149149                  ! horizontal diffusive trends 
    150                   zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     150                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    151151                  ztra = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) 
    152152                  ! add it to the general tracer trends 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r2715 r2905  
    172172            DO jj = 1 , jpjm1 
    173173               DO ji = 1, fs_jpim1   ! vector opt. 
    174                   zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) 
    175                   zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) 
     174                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e1ur(ji,jj) * fse3u_n(ji,jj,jk) 
     175                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e2vr(ji,jj) * fse3v_n(ji,jj,jk) 
    176176                  ! 
    177177                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   & 
    … …  
    197197            DO jj = 2 , jpjm1 
    198198               DO ji = fs_2, fs_jpim1   ! vector opt. 
    199                   zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     199                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    200200                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
    201201                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    … …  
    282282            DO jj = 2, jpjm1 
    283283               DO ji = fs_2, fs_jpim1   ! vector opt. 
    284                   zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     284                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    285285                  ztra = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr 
    286286                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r2715 r2905  
    3030 
    3131   PUBLIC   tra_ldf_lap   ! routine called by step.F90 
    32  
    33    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::   e1ur, e2vr   ! scale factor coefficients 
    3432 
    3533   !! * Substitutions 
    … …  
    8179         IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype 
    8280         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    83          ! 
    84          ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr ) 
    85          IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    86          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' ) 
    87          ! 
    88          e1ur(:,:) = e2u(:,:) / e1u(:,:) 
    89          e2vr(:,:) = e1v(:,:) / e2v(:,:) 
    9081      ENDIF 
    9182 
    … …  
    9990            DO jj = 1, jpjm1 
    10091               DO ji = 1, fs_jpim1   ! vector opt. 
    101                   zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u(ji,jj,jk) 
    102                   zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v(ji,jj,jk) 
     92                  zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u_n(ji,jj,jk) 
     93                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v_n(ji,jj,jk) 
    10394                  ztu(ji,jj,jk) = zabe1 * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
    10495                  ztv(ji,jj,jk) = zabe2 * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
    … …  
    112103                     ikv = mbkv(ji,jj) 
    113104                     IF( iku == jk ) THEN 
    114                         zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku) 
     105                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u_n(ji,jj,iku) 
    115106                        ztu(ji,jj,jk) = zabe1 * pgu(ji,jj,jn) 
    116107                     ENDIF 
    117108                     IF( ikv == jk ) THEN 
    118                         zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v(ji,jj,ikv) 
     109                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v_n(ji,jj,ikv) 
    119110                        ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 
    120111                     ENDIF 
    … …  
    128119            DO jj = 2, jpjm1 
    129120               DO ji = fs_2, fs_jpim1   ! vector opt. 
    130                   zbtr = 1._wp / ( e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     121                  zbtr = 1._wp / ( e1t(ji,jj) *e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    131122                  ! horizontal diffusive trends added to the general tracer trends 
    132123                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r2715 r2905  
    3636   !! ------------                                      !  fields  ! fields ! trends ! 
    3737   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m] 
    38    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m] 
    39    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m] 
    40    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::            sshf_n          !: sea surface height at f-point [m] 
    4138   ! 
    4239   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   spgu, spgv               !: horizontal surface pressure gradient 
    … …  
    7471         &     rhop(jpi,jpj,jpk) ,                                         & 
    7572         &     sshb  (jpi,jpj)   , sshn  (jpi,jpj) , ssha  (jpi,jpj) ,     & 
    76          &     sshu_b(jpi,jpj)   , sshu_n(jpi,jpj) , sshu_a(jpi,jpj) ,     & 
    77          &     sshv_b(jpi,jpj)   , sshv_n(jpi,jpj) , sshv_a(jpi,jpj) ,     & 
    78          &                         sshf_n(jpi,jpj) ,                       & 
    7973         &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       & 
    8074         &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     & 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90

    r2715 r2905  
    102102 
    103103      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    104       !  Ocean dynamics : ssh, wn, hdiv, rot                                 ! 
    105       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    106                          CALL ssh_wzv( kstp )         ! after ssh & vertical velocity 
    107  
     104      !  Ocean dynamics : hdiv, rot, ssh, e3, wn 
     105      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     106                         CALL ssh_nxt       ( kstp )  ! after ssh 
     107      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
     108                         CALL wzv           ( kstp )  ! now cross-level velocity 
    108109      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    109110      ! Ocean physics update                (ua, va, ta, sa used as workspace) 
    … …  
    232233                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step 
    233234 
    234                                CALL ssh_nxt( kstp )         ! sea surface height at next time step 
     235                               CALL ssh_swp( kstp )         ! swap of sea surface height 
     236      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
    235237 
    236238      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics 
    237       IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
     239      IF( lk_diaobs        )   CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
    238240 
    239241      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 

  • branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r2528 r2905  
    5757 
    5858   USE sshwzv           ! vertical velocity and ssh        (ssh_wzv routine) 
     59   USE domvvl           ! variable vertical scale factors  (dom_vvl_sf_nxt routine) 
     60   !                                                       (dom_vvl_sf_swp routine) 
    5961 
    6062   USE ldfslp           ! iso-neutral slopes               (ldf_slp routine) 
Note: See TracChangeset for help on using the changeset viewer.