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 1438 – NEMO

Changeset 1438


Ignore:
Timestamp:
2009-05-11T16:34:47+02:00 (15 years ago)
Author:
rblod
Message:

Merge VVL branch with the trunk (act II), see ticket #429

Location:
trunk/NEMO
Files:
24 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMO/NST_SRC/agrif_opa_update.F90

    r1300 r1438  
    8080 
    8181 
    82       CALL Agrif_ChildGrid_To_ParentGrid() 
    83       CALL recompute_diags( kt ) 
    84       CALL Agrif_ParentGrid_To_ChildGrid() 
     82!Done in step 
     83!      CALL Agrif_ChildGrid_To_ParentGrid() 
     84!      CALL recompute_diags( kt ) 
     85!      CALL Agrif_ParentGrid_To_ChildGrid() 
    8586 
    8687#endif 

  • trunk/NEMO/OPA_SRC/DOM/dom_oce.F90

    r1241 r1438  
    114114#if defined key_vvl 
    115115   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .TRUE.    !: variable grid flag 
    116  
    117116   !! All coordinates 
    118117   !! --------------- 
    119    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !: 
    120       gdep3w_1        ,  &  !: depth of T-points (sum of e3w) (m) 
    121       gdept_1, gdepw_1,  &  !: analytical depth at T-W  points (m) 
    122       e3v_1  , e3f_1  ,  &  !: analytical vertical scale factors at  V--F 
    123       e3t_1  , e3u_1  ,  &  !:                                       T--U  points (m) 
    124       e3vw_1          ,  &  !: analytical vertical scale factors at  VW-- 
    125       e3w_1  , e3uw_1       !:                                       W--UW  points (m) 
     118   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   gdep3w_1           !: depth of T-points (sum of e3w) (m) 
     119   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   gdept_1, gdepw_1   !: analytical depth at T-W  points (m) 
     120   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3v_1  , e3f_1     !: analytical vertical scale factors at  V--F 
     121   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_1  , e3u_1     !:                                       T--U  points (m) 
     122   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3vw_1             !: analytical vertical scale factors at  VW-- 
     123   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3w_1  , e3uw_1    !:                                       W--UW  points (m) 
    126124#else 
    127125   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag 
    … …  
    129127   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !: 
    130128      hur, hvr,          &  !: inverse of u and v-points ocean depth (1/m) 
    131       hu , hv               !: depth at u- and v-points (meters) 
     129      hu , hv,           &  !: depth at u- and v-points (meters) 
     130      hu_0 , hv_0           !: refernce depth at u- and v-points (meters) 
    132131 
    133132   !! z-coordinate with full steps (also used in the other cases as reference z-coordinate) 

  • trunk/NEMO/OPA_SRC/DOM/domain.F90

    r1335 r1438  
    44   !! Ocean initialization : domain initialization 
    55   !!============================================================================== 
    6  
     6   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code 
     7   !!                 !  1991-11  (G. Madec) 
     8   !!                 !  1992-01  (M. Imbard) insert time step initialization 
     9   !!                 !  1996-06  (G. Madec) generalized vertical coordinate  
     10   !!                 !  1997-02  (G. Madec) creation of domwri.F 
     11   !!                 !  2001-05  (E.Durand - G. Madec) insert closed sea 
     12   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module 
     13   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
     14   !!---------------------------------------------------------------------- 
     15    
    716   !!---------------------------------------------------------------------- 
    817   !!   dom_init       : initialize the space and time domain 
    … …  
    1019   !!   dom_ctl        : control print for the ocean domain 
    1120   !!---------------------------------------------------------------------- 
    12    !! * Modules used 
    1321   USE oce             !  
    1422   USE dom_oce         ! ocean space and time domain 
    … …  
    3038   PRIVATE 
    3139 
    32    !! * Routine accessibility 
    33    PUBLIC dom_init       ! called by opa.F90 
     40   PUBLIC   dom_init   ! called by opa.F90 
    3441 
    3542   !! * Substitutions 
    3643#  include "domzgr_substitute.h90" 
    37    !!---------------------------------------------------------------------- 
    38    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     44   !!------------------------------------------------------------------------- 
     45   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    3946   !! $Id$ 
    40    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    41    !!---------------------------------------------------------------------- 
     47   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     48   !!------------------------------------------------------------------------- 
    4249 
    4350CONTAINS 
    … …  
    5865      !!      - dom_stp: defined the model time step 
    5966      !!      - dom_wri: create the meshmask file if nmsh=1 
    60       !! 
    61       !! History : 
    62       !!        !  90-10  (C. Levy - G. Madec)  Original code 
    63       !!        !  91-11  (G. Madec) 
    64       !!        !  92-01  (M. Imbard) insert time step initialization 
    65       !!        !  96-06  (G. Madec) generalized vertical coordinate  
    66       !!        !  97-02  (G. Madec) creation of domwri.F 
    67       !!        !  01-05  (E.Durand - G. Madec) insert closed sea 
    68       !!   8.5  !  02-08  (G. Madec)  F90: Free form and module 
    69       !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
    70       !!---------------------------------------------------------------------- 
    71       !! * Local declarations 
     67      !!---------------------------------------------------------------------- 
    7268      INTEGER ::   jk                ! dummy loop argument 
    7369      INTEGER ::   iconf = 0         ! temporary integers 
    … …  
    9086      CALL dom_msk                        ! Masks 
    9187 
    92       IF( lk_vvl )   CALL dom_vvl_ini     ! Vertical variable mesh 
     88      IF( lk_vvl )   CALL dom_vvl         ! Vertical variable mesh 
    9389 
    9490      ! Local depth or Inverse of the local depth of the water column at u- and v-points 
    9591      ! ------------------------------ 
    9692      ! Ocean depth at U- and V-points 
    97       hu(:,:) = 0. 
    98       hv(:,:) = 0. 
    99  
     93      hu(:,:) = 0.e0 
     94      hv(:,:) = 0.e0 
    10095      DO jk = 1, jpk 
    10196         hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    … …  
    105100      hur(:,:) = fse3u(:,:,1)             ! Lower bound : thickness of the first model level 
    106101      hvr(:,:) = fse3v(:,:,1) 
    107  
    108102      DO jk = 2, jpk                      ! Sum of the vertical scale factors 
    109103         hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    110104         hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
    111105      END DO 
    112  
    113106      ! Compute and mask the inverse of the local depth 
    114107      hur(:,:) = 1. / hur(:,:) * umask(:,:,1) 
    115108      hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1) 
    116109 
    117  
    118110      CALL dom_stp                        ! Time step 
    119111 
    … …  
    121113 
    122114      IF( .NOT.ln_rstart )   CALL dom_ctl    ! Domain control 
    123  
     115      ! 
    124116   END SUBROUTINE dom_init 
    125117 
    … …  
    134126      !!              - namdom namelist 
    135127      !!              - namcla namelist 
    136       !! 
    137       !! History : 
    138       !!   9.0  !  03-08  (G. Madec)  Original code 
    139       !!---------------------------------------------------------------------- 
    140       !! * Modules used 
     128      !!---------------------------------------------------------------------- 
    141129      USE ioipsl 
    142130      NAMELIST/namrun/ no    , cexper, cn_ocerst_in, cn_ocerst_out, ln_rstart, nrstdt,   & 
    … …  
    156144      ENDIF 
    157145 
    158       ! Namelist namrun : parameters of the run 
    159       REWIND( numnam ) 
     146      REWIND( numnam )              ! Namelist namrun : parameters of the run 
    160147      READ  ( numnam, namrun ) 
    161  
    162148      IF(lwp) THEN 
    163149         WRITE(numout,*) '        Namelist namrun' 
    … …  
    228214      ENDIF 
    229215 
    230       ! Namelist namdom : space/time domain (bathymetry, mesh, timestep) 
    231       REWIND( numnam ) 
     216      REWIND( numnam )              ! Namelist namdom : space/time domain (bathymetry, mesh, timestep) 
    232217      READ  ( numnam, namdom ) 
    233218 
    … …  
    252237      ENDIF 
    253238 
    254       ! Default values 
    255239      n_cla = 0 
    256  
    257       ! Namelist cross land advection 
    258       REWIND( numnam ) 
     240      REWIND( numnam )              ! Namelist cross land advection 
    259241      READ  ( numnam, namcla ) 
    260242      IF(lwp) THEN 
    … …  
    264246      ENDIF 
    265247 
    266       IF( nbit_cmp == 1 .AND. n_cla /= 0 ) THEN 
    267          CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require n_cla = 0' ) 
    268       END IF 
    269  
     248      IF( nbit_cmp == 1 .AND. n_cla /= 0 )   CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require n_cla = 0' ) 
     249      ! 
    270250   END SUBROUTINE dom_nam 
    271251 
    … …  
    278258      !! 
    279259      !! ** Method  :   compute and print extrema of masked scale factors 
    280       !! 
    281       !! History : 
    282       !!   8.5  !  02-08  (G. Madec)    Original code 
    283       !!---------------------------------------------------------------------- 
    284       !! * Local declarations 
     260      !!---------------------------------------------------------------------- 
    285261      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2 
    286262      INTEGER, DIMENSION(2) ::   iloc      !  
    287263      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max 
    288264      !!---------------------------------------------------------------------- 
    289  
    290       ! Extrema of the scale factors 
    291265 
    292266      IF(lwp)WRITE(numout,*) 
    … …  
    325299         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 
    326300      ENDIF 
    327  
     301      ! 
    328302   END SUBROUTINE dom_ctl 
    329303 

  • trunk/NEMO/OPA_SRC/DOM/domvvl.F90

    r1146 r1438  
    44   !! Ocean :  
    55   !!====================================================================== 
    6    !! History :   9.0  !  06-06  (B. Levier, L. Marie)  original code 
    7    !!              "   !  07-07  (D. Storkey) Bug fixes and code for BDY option. 
     6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code 
     7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate 
    88   !!---------------------------------------------------------------------- 
    9  
     9#if defined key_vvl 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   'key_vvl'                              variable volume 
    1212   !!---------------------------------------------------------------------- 
     13   !!   dom_vvl     : defined coefficients to distribute ssh on each layers 
    1314   !!---------------------------------------------------------------------- 
    14    !!   dom_vvl         : defined scale factors & depths at each time step 
    15    !!   dom_vvl_ini     : defined coefficients to distribute ssh on each layers 
    16    !!   dom_vvl_ssh     : defined the ocean sea level at each time step 
    17    !!---------------------------------------------------------------------- 
    18    !! * Modules used 
    1915   USE oce             ! ocean dynamics and tracers 
    2016   USE dom_oce         ! ocean space and time domain 
    2117   USE sbc_oce         ! surface boundary condition: ocean 
    22    USE dynspg_oce      ! surface pressure gradient variables 
    2318   USE phycst          ! physical constants 
    2419   USE in_out_manager  ! I/O manager 
    2520   USE lib_mpp         ! distributed memory computing library 
    2621   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    27    USE bdy_oce         ! unstructured open boundary conditions 
    2822 
    2923   IMPLICIT NONE 
    3024   PRIVATE 
    3125 
    32    !! * Routine accessibility 
    33    PUBLIC dom_vvl_ini    ! called by dom_init.F90 
    34    PUBLIC dom_vvl_ssh    ! called by trazdf.F90 
    35    PUBLIC dom_vvl        ! called by istate.F90 and step.F90 
    36    PUBLIC dom_vvl_sf_ini !  
    37    PUBLIC dom_vvl_sf     !  
     26   PUBLIC dom_vvl        ! called by domain.F90 
    3827 
    39    !! * Module variables 
    40    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !: 
    41       mut, muu, muv, muf                            !: 
     28   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ee_t, ee_u, ee_v, ee_f   !: ??? 
     29 
     30   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   mut, muu, muv, muf   !: ???  
    4231 
    4332   REAL(wp), DIMENSION(jpk) ::   r2dt               ! vertical profile time-step, = 2 rdttra  
    … …  
    4837#  include "vectopt_loop_substitute.h90" 
    4938   !!---------------------------------------------------------------------- 
    50    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     39   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    5140   !! $Id$ 
    5241   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    5544CONTAINS        
    5645 
    57 #if defined key_vvl 
    58  
    59    SUBROUTINE dom_vvl_ini 
     46   SUBROUTINE dom_vvl 
    6047      !!---------------------------------------------------------------------- 
    61       !!                ***  ROUTINE dom_vvl_ini  *** 
     48      !!                ***  ROUTINE dom_vvl  *** 
    6249      !!                    
    6350      !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread 
    … …  
    6552      !! 
    6653      !!---------------------------------------------------------------------- 
    67       INTEGER  ::   ji, jj, jk, zpk 
    68       REAL(wp), DIMENSION(jpi,jpj) ::   zmut, zmuu, zmuv, zmuf   ! 2D workspace 
     54      INTEGER  ::   ji, jj, jk 
     55      REAL(wp) ::   zcoefu, zcoefv, zcoeff 
    6956      !!---------------------------------------------------------------------- 
    7057 
    7158      IF(lwp)   THEN 
    7259         WRITE(numout,*) 
    73          WRITE(numout,*) 'dom_vvl_ini : Variable volume activated' 
    74          WRITE(numout,*) '~~~~~~~~~~~   compute coef. used to spread ssh over each layers' 
     60         WRITE(numout,*) 'dom_vvl : Variable volume activated' 
     61         WRITE(numout,*) '~~~~~~~~  compute coef. used to spread ssh over each layers' 
    7562      ENDIF 
    76  
    77       IF( ln_zps )  CALL ctl_stop( 'dom_vvl_ini : option ln_zps is incompatible with variable volume option key_vvl') 
    7863 
    7964#if defined key_zco  ||  defined key_dynspg_rl 
    … …  
    8166#endif 
    8267 
    83       fsvdept (:,:,:) = gdept (:,:,:) 
    84       fsvdepw (:,:,:) = gdepw (:,:,:) 
    85       fsvde3w (:,:,:) = gdep3w(:,:,:) 
    86       fsve3t (:,:,:) = e3t   (:,:,:) 
    87       fsve3u (:,:,:) = e3u   (:,:,:) 
    88       fsve3v (:,:,:) = e3v   (:,:,:) 
    89       fsve3f (:,:,:) = e3f   (:,:,:) 
    90       fsve3w (:,:,:) = e3w   (:,:,:) 
    91       fsve3uw (:,:,:) = e3uw  (:,:,:) 
    92       fsve3vw (:,:,:) = e3vw  (:,:,:) 
     68      fsdept(:,:,:) = gdept (:,:,:) 
     69      fsdepw(:,:,:) = gdepw (:,:,:) 
     70      fsde3w(:,:,:) = gdep3w(:,:,:) 
     71      fse3t (:,:,:) = e3t   (:,:,:) 
     72      fse3u (:,:,:) = e3u   (:,:,:) 
     73      fse3v (:,:,:) = e3v   (:,:,:) 
     74      fse3f (:,:,:) = e3f   (:,:,:) 
     75      fse3w (:,:,:) = e3w   (:,:,:) 
     76      fse3uw(:,:,:) = e3uw  (:,:,:) 
     77      fse3vw(:,:,:) = e3vw  (:,:,:) 
    9378 
    9479      ! mu computation 
    95       zmut(:,:)   = 0.e0 
    96       zmuu(:,:)   = 0.e0 
    97       zmuv(:,:)   = 0.e0 
    98       zmuf(:,:)   = 0.e0 
    99       mut (:,:,:) = 0.e0 
    100       muu (:,:,:) = 0.e0 
    101       muv (:,:,:) = 0.e0 
    102       muf (:,:,:) = 0.e0 
     80      ! -------------- 
     81      ! define ee_t, u, v and f as in sigma coordinate (ee_t = 1/ht, ...) 
     82      ee_t(:,:) = fse3t_0(:,:,1)        ! Lower bound : thickness of the first model level 
     83      ee_u(:,:) = fse3u_0(:,:,1) 
     84      ee_v(:,:) = fse3v_0(:,:,1) 
     85      ee_f(:,:) = fse3f_0(:,:,1) 
     86      DO jk = 2, jpkm1                   ! Sum of the masked vertical scale factors 
     87         ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 
     88         ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
     89         ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
     90         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask 
     91            ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk) 
     92         END DO 
     93      END DO   
     94      !                                  ! Compute and mask the inverse of the local depth at T, U, V and F points 
     95      ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1) 
     96      ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1) 
     97      ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1) 
     98      DO jj = 1, jpjm1                         ! f-point case fmask cannot be used  
     99         ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 
     100      END DO 
     101      CALL lbc_lnk( ee_f, 'F', 1. )       ! lateral boundary condition on ee_f 
     102      ! 
     103      DO jk = 1, jpk 
     104         mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk)   ! at T levels 
     105         muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk)   ! at T levels 
     106         muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk)   ! at T levels 
     107      END DO 
     108      DO jk = 1, jpk 
     109         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask 
     110               muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels 
     111         END DO 
     112         muf(:,jpj,jk) = 0.e0 
     113      END DO 
     114      CALL lbc_lnk( muf, 'F', 1. )       ! lateral boundary condition on ee_f 
    103115 
    104       DO jj = 1, jpj 
    105          DO ji = 1, jpi 
    106             zpk = mbathy(ji,jj) - 1 
    107             DO jk = 1, zpk 
    108                zmut(ji,jj) = zmut(ji,jj) + fsve3t(ji,jj,jk) * SUM( fsve3t(ji,jj,jk:zpk) ) 
    109                zmuu(ji,jj) = zmuu(ji,jj) + fsve3u(ji,jj,jk) * SUM( fsve3u(ji,jj,jk:zpk) ) 
    110                zmuv(ji,jj) = zmuv(ji,jj) + fsve3v(ji,jj,jk) * SUM( fsve3v(ji,jj,jk:zpk) ) 
    111                zmuf(ji,jj) = zmuf(ji,jj) + fsve3f(ji,jj,jk) * SUM( fsve3f(ji,jj,jk:zpk) ) 
    112             END DO 
     116 
     117      ! Reference ocean depth at U- and V-points 
     118      hu_0(:,:) = 0.e0     
     119      hv_0(:,:) = 0.e0 
     120      DO jk = 1, jpk 
     121         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 
     122         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 
     123      END DO 
     124 
     125      ! before and now Sea Surface Height at u-, v-, f-points 
     126      DO jj = 1, jpjm1 
     127         DO ji = 1, jpim1 
     128            zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     129            zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     130            zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
     131            sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
     132               &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )    
     133            sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
     134               &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )    
     135            sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 & 
     136               &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) )                
     137            sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
     138               &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )    
     139            sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &  
     140               &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )      
     141            sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 & 
     142               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) )                
    113143         END DO 
    114144      END DO 
    115  
    116       DO jj = 1, jpj 
    117          DO ji = 1, jpi 
    118             zpk = mbathy(ji,jj) - 1 
    119             DO jk = 1, zpk 
    120 #if defined key_sigma_vvl 
    121                mut(ji,jj,jk) = 1./SUM( fsve3t(ji,jj,1:zpk) )  
    122                muu(ji,jj,jk) = 1./SUM( fsve3u(ji,jj,1:zpk) )  
    123                muv(ji,jj,jk) = 1./SUM( fsve3v(ji,jj,1:zpk) )  
    124                muf(ji,jj,jk) = 1./SUM( fsve3f(ji,jj,1:zpk) )  
    125 #else 
    126                mut(ji,jj,jk) = SUM( fsve3t(ji,jj,jk:zpk) ) / zmut(ji,jj) 
    127                muu(ji,jj,jk) = SUM( fsve3u(ji,jj,jk:zpk) ) / zmuu(ji,jj) 
    128                muv(ji,jj,jk) = SUM( fsve3v(ji,jj,jk:zpk) ) / zmuv(ji,jj) 
    129                muf(ji,jj,jk) = SUM( fsve3f(ji,jj,jk:zpk) ) / zmuf(ji,jj) 
    130 #endif 
    131             END DO 
    132          END DO 
    133       END DO 
    134  
    135  
    136    END SUBROUTINE dom_vvl_ini 
    137  
    138  
    139  
    140    SUBROUTINE dom_vvl 
    141       !!---------------------------------------------------------------------- 
    142       !!                ***  ROUTINE dom_vvl  *** 
    143       !!                    
    144       !! ** Purpose :  compute ssh at U-V-F points, T-W scale factors and local 
    145       !!               depths at each time step. 
    146       !!---------------------------------------------------------------------- 
    147       !! * Local declarations 
    148       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    149       REAL(wp), DIMENSION(jpi,jpj) :: zsshf    ! 2D workspace 
    150       !!---------------------------------------------------------------------- 
    151  
    152       ! Sea Surface Height at u-v-fpoints 
    153       DO jj = 1, jpjm1 
    154          DO ji = 1,jpim1 
    155             sshu(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  & 
    156                &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    157                &            + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
    158             ! 
    159             sshv(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  & 
    160                &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
    161                &            + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
    162             ! 
    163             zsshf(ji,jj) = 0.25 * fmask(ji,jj,1)                 & 
    164                &           * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
    165                &             + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 
    166          END DO 
    167       END DO 
    168  
    169145      ! Boundaries conditions 
    170       CALL lbc_lnk( sshu,  'U', 1. ) 
    171       CALL lbc_lnk( sshv,  'V', 1. ) 
    172       CALL lbc_lnk( zsshf, 'F', 1. ) 
    173  
    174       ! Scale factors at T levels 
    175       DO jk = 1, jpkm1 
    176          fse3t(:,:,jk) = fsve3t(:,:,jk) * ( 1 + sshn(:,:)  * mut(:,:,jk) ) 
    177          fse3u(:,:,jk) = fsve3u(:,:,jk) * ( 1 + sshu(:,:)  * muu(:,:,jk) ) 
    178          fse3v(:,:,jk) = fsve3v(:,:,jk) * ( 1 + sshv(:,:)  * muv(:,:,jk) ) 
    179          fse3f(:,:,jk) = fsve3f(:,:,jk) * ( 1 + zsshf(:,:) * muf(:,:,jk) ) 
    180       END DO 
    181  
    182       ! Scale factors at W levels 
    183       fse3w (:,:,1) = fse3t(:,:,1) 
    184       fse3uw(:,:,1) = fse3u(:,:,1) 
    185       fse3vw(:,:,1) = fse3v(:,:,1) 
    186       DO jk = 2, jpk 
    187          fse3w (:,:,jk) = 0.5 * ( fse3t(:,:,jk-1) + fse3t(:,:,jk) ) 
    188          fse3uw(:,:,jk) = 0.5 * ( fse3u(:,:,jk-1) + fse3u(:,:,jk) ) 
    189          fse3vw(:,:,jk) = 0.5 * ( fse3v(:,:,jk-1) + fse3v(:,:,jk) ) 
    190       END DO 
    191  
    192       ! T and W points depth 
    193       fsdept(:,:,1) = 0.5 * fse3w(:,:,1) 
    194       fsdepw(:,:,1) = 0.e0 
    195       fsde3w(:,:,1) = fsdept(:,:,1) - sshn(:,:) 
    196       DO jk = 2, jpk 
    197          fsdept(:,:,jk) = fsdept(:,:,jk-1) + fse3w(:,:,jk) 
    198          fsdepw(:,:,jk) = fsdepw(:,:,jk-1) + fse3t(:,:,jk-1) 
    199          fsde3w(:,:,jk) = fsdept(:,:,jk  ) - sshn(:,:) 
    200       END DO 
    201  
    202       IF( MINVAL(fse3t(:,:,:)) < 0.0 ) THEN 
    203          CALL ctl_stop('E R R O R : dom_vvl','Level thickness fse3t less than zero.') 
    204          nstop = nstop+1 
    205       ENDIF 
    206  
    207       ! Local depth or Inverse of the local depth of the water column at u- and v-points 
    208       ! ------------------------------ 
    209  
    210       ! Ocean depth at U- and V-points 
    211       hu(:,:) = 0.e0    ;    hv(:,:) = 0.e0 
    212  
    213       DO jk = 1, jpk 
    214          hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    215          hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
    216       END DO 
    217  
    218  
    219       ! Inverse of the local depth 
    220       hur(:,:) = fse3u(:,:,1)             ! Lower bound : thickness of the first model level 
    221       hvr(:,:) = fse3v(:,:,1) 
    222        
    223       DO jk = 2, jpk                      ! Sum of the vertical scale factors 
    224          hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 
    225          hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 
    226       END DO 
    227  
    228       ! Compute and mask the inverse of the local depth 
    229       hur(:,:) = 1. / hur(:,:) * umask(:,:,1) 
    230       hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1) 
    231  
     146      CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
     147      CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
     148      CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. ) 
     149      ! 
    232150   END SUBROUTINE dom_vvl 
    233151 
    234  
    235  
    236    SUBROUTINE dom_vvl_ssh( kt )  
    237       !!---------------------------------------------------------------------- 
    238       !!                ***  ROUTINE dom_vvl_ssh  *** 
    239       !!                    
    240       !! ** Purpose :  compute the ssha field used in tra_zdf, dynnxt, tranxt  
    241       !!               and dynspg routines 
    242       !!---------------------------------------------------------------------- 
    243       !! * Arguments 
    244       INTEGER, INTENT( in )  ::    kt                         ! time step 
    245  
    246       !! * Local declarations 
    247       INTEGER  ::   ji, jj, jk                                ! dummy loop indices 
    248       REAL(wp) ::   z2dt, zraur                               ! temporary scalars 
    249       REAL(wp), DIMENSION(jpi,jpj) ::   zun, zvn, zhdiv       ! 2D workspace 
    250       !!---------------------------------------------------------------------- 
    251  
    252       !! Sea surface height after (ssha) in variable volume case 
    253       !! --------------------------====-----===============----- 
    254       !! ssha is needed for the tracer time-stepping (trazdf_imp or 
    255       !! tra_nxt), for the ssh time-stepping (dynspg_*) and for the 
    256       !! dynamics time-stepping (dynspg_flt or dyn_nxt, and wzv). 
    257       !! un, vn (or un_b and vn_b) and emp are needed, so it must be 
    258       !! done after the call of oce_sbc. 
    259  
    260       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000 
    261          r2dt(:) =  rdttra(:)                       ! = rdtra (restarting with Euler time stepping) 
    262       ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1 
    263          r2dt(:) = 2. * rdttra(:)                   ! = 2 rdttra (leapfrog) 
    264       ENDIF 
    265  
    266       z2dt  = r2dt(1) 
    267       zraur = 1. / rauw 
    268  
    269       ! Vertically integrated quantities 
    270       ! -------------------------------- 
    271       IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN 
    272          zun(:,:) = 0.e0    ;    zvn(:,:) = 0.e0 
    273          ! 
    274          DO jk = 1, jpkm1 
    275             !                                               ! Vertically integrated transports (now) 
    276             zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    277             zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    278          END DO 
    279       ELSE ! lk_dynspg_ts is true 
    280          zun (:,:) = un_b(:,:)    ;    zvn (:,:) = vn_b(:,:) 
    281       ENDIF 
    282  
    283       ! Horizontal divergence of barotropic transports 
    284       !-------------------------------------------------- 
    285       DO jj = 2, jpjm1 
    286          DO ji = 2, jpim1   ! vector opt. 
    287             zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )    & 
    288                &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )    & 
    289                &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )    & 
    290                &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )  & 
    291                &           / ( e1t(ji,jj) * e2t(ji,jj) ) 
    292          END DO 
    293       END DO 
    294  
    295 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 
    296       ! open boundaries (div must be zero behind the open boundary) 
    297       !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    298       IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east 
    299       IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west 
    300       IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north 
    301       IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south 
    302 #endif 
    303  
    304 #if defined key_bdy 
    305       zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) 
    306 #endif 
    307  
    308       CALL lbc_lnk( zhdiv, 'T', 1. ) 
    309  
    310       ! Sea surface elevation time stepping 
    311       ! ----------------------------------- 
    312       IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN 
    313          ssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 
    314       ELSE ! lk_dynspg_ts is true 
    315          ssha(:,:) = (  sshb_b(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    316       ENDIF 
    317  
    318  
    319    END SUBROUTINE dom_vvl_ssh 
    320  
    321    SUBROUTINE  dom_vvl_sf( zssh, gridp, sfe3 ) 
    322       !!---------------------------------------------------------------------- 
    323       !!                ***  ROUTINE dom_vvl_sf  *** 
    324       !!                    
    325       !! ** Purpose :  compute vertical scale factor at each time step 
    326       !!---------------------------------------------------------------------- 
    327       !! * Arguments 
    328       CHARACTER(LEN=1)                , INTENT( in ) :: gridp   ! grid point type 
    329       REAL(wp), DIMENSION(jpi,jpj)    , INTENT( in ) :: zssh    ! 2D workspace 
    330       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: sfe3    ! 3D workspace 
    331  
    332       !! * Local declarations 
    333       INTEGER  ::   jk                                      ! dummy loop indices 
    334       !!---------------------------------------------------------------------- 
    335  
    336       SELECT CASE (gridp) 
    337  
    338       CASE ('U') 
    339          ! 
    340          DO jk = 1, jpk 
    341             sfe3(:,:,jk)  = fsve3u(:,:,jk) * ( 1 + zssh(:,:) * muu(:,:,jk) ) 
    342          ENDDO 
    343  
    344       CASE ('V') 
    345          ! 
    346          DO jk = 1, jpk 
    347             sfe3(:,:,jk)  = fsve3v(:,:,jk) * ( 1 + zssh(:,:) * muv(:,:,jk) ) 
    348          ENDDO 
    349  
    350       CASE ('T') 
    351          ! 
    352          DO jk = 1, jpk 
    353             sfe3(:,:,jk)  = fsve3t(:,:,jk) * ( 1 + zssh(:,:) * mut(:,:,jk) ) 
    354          ENDDO 
    355  
    356       END SELECT 
    357  
    358    END SUBROUTINE dom_vvl_sf 
    359  
    360    SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini ) 
    361       !!---------------------------------------------------------------------- 
    362       !!                ***  ROUTINE dom_vvl_sf_ini  *** 
    363       !!                    
    364       !! ** Purpose :  affect the appropriate vertical scale factor. It is done 
    365       !!               to avoid compiling error when using key_zco cpp key 
    366       !!---------------------------------------------------------------------- 
    367       !! * Arguments 
    368       CHARACTER(LEN=1)                , INTENT( in ) :: gridp      ! grid point type 
    369       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini    ! 3D workspace 
    370       !!---------------------------------------------------------------------- 
    371  
    372       SELECT CASE (gridp) 
    373  
    374       CASE ('U') 
    375          ! 
    376          sfe3ini(:,:,:)  = fse3u(:,:,:) 
    377  
    378       CASE ('V') 
    379          ! 
    380          sfe3ini(:,:,:)  = fse3v(:,:,:) 
    381  
    382       CASE ('T') 
    383          ! 
    384          sfe3ini(:,:,:)  = fse3t(:,:,:) 
    385  
    386       END SELECT 
    387  
    388    END SUBROUTINE dom_vvl_sf_ini 
    389152#else 
    390153   !!---------------------------------------------------------------------- 
    391154   !!   Default option :                                      Empty routine 
    392155   !!---------------------------------------------------------------------- 
    393    SUBROUTINE dom_vvl_ini 
    394    END SUBROUTINE dom_vvl_ini 
     156CONTAINS 
    395157   SUBROUTINE dom_vvl 
    396158   END SUBROUTINE dom_vvl 
    397    SUBROUTINE dom_vvl_ssh( kt )  
    398       !! * Arguments 
    399       INTEGER, INTENT( in )  ::    kt        ! time step 
    400       WRITE(*,*) 'dom_vvl_ssh: You should not have seen this print! error?', kt 
    401    END SUBROUTINE dom_vvl_ssh 
    402    SUBROUTINE dom_vvl_sf( zssh, gridp, sfe3 )  
    403       !! * Arguments 
    404       CHARACTER(LEN=1)                , INTENT( in ) :: gridp   ! grid point type 
    405       REAL(wp), DIMENSION(jpi,jpj)    , INTENT( in ) :: zssh    ! 2D workspace 
    406       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3    ! 3D workspace 
    407       sfe3(:,:,:) = 0.e0 
    408       WRITE(*,*) 'sfe3: You should not have seen this print! error?', gridp 
    409       WRITE(*,*) 'sfe3: You should not have seen this print! error?', zssh(1,1) 
    410    END SUBROUTINE dom_vvl_sf 
    411    SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini )  
    412       !! * Arguments 
    413       CHARACTER(LEN=1)                , INTENT( in ) :: gridp    ! grid point type 
    414       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini  ! 3D workspace 
    415       sfe3ini(:,:,:) = 0.e0 
    416       WRITE(*,*) 'sfe3ini: You should not have seen this print! error?', gridp 
    417    END SUBROUTINE dom_vvl_sf_ini 
    418159#endif 
    419160 

  • trunk/NEMO/OPA_SRC/DOM/domzgr_substitute.h90

    r1156 r1438  
    55   !!      factors depending on the vertical coord. used, using CPP macro. 
    66   !!---------------------------------------------------------------------- 
     7   !! History :  1.0  !  2005-10  (A. Beckmann, G. Madec) generalisation to all coord. 
     8   !!            3.1  !  2009-02  (G. Madec, M. Leclair)  pure z* coordinate 
    79   !!---------------------------------------------------------------------- 
    8    !!  OPA 9.0 , LOCEAN-IPSL (2005)  
     10#if defined key_zco 
     11! reference for pure z-coordinate (1D - no i,j and time dependency) 
     12#   define  fsdept_0(i,j,k)  gdept_0(k) 
     13#   define  fsdepw_0(i,j,k)  gdepw_0(k) 
     14#   define  fsde3w_0(i,j,k)  gdepw_0(k) 
     15#   define  fse3t_0(i,j,k)   e3t_0(k) 
     16#   define  fse3u_0(i,j,k)   e3t_0(k) 
     17#   define  fse3v_0(i,j,k)   e3t_0(k) 
     18#   define  fse3f_0(i,j,k)   e3t_0(k) 
     19#   define  fse3w_0(i,j,k)   e3w_0(k) 
     20#   define  fse3uw_0(i,j,k)  e3w_0(k) 
     21#   define  fse3vw_0(i,j,k)  e3w_0(k) 
     22#else 
     23! reference for s- or zps-coordinate (3D no time dependency) 
     24#   define  fsdept_0(i,j,k)  gdept(i,j,k) 
     25#   define  fsdepw_0(i,j,k)  gdepw(i,j,k) 
     26#   define  fsde3w_0(i,j,k)  gdep3w(i,j,k) 
     27#   define  fse3t_0(i,j,k)   e3t(i,j,k) 
     28#   define  fse3u_0(i,j,k)   e3u(i,j,k) 
     29#   define  fse3v_0(i,j,k)   e3v(i,j,k) 
     30#   define  fse3f_0(i,j,k)   e3f(i,j,k) 
     31#   define  fse3w_0(i,j,k)   e3w(i,j,k) 
     32#   define  fse3uw_0(i,j,k)  e3uw(i,j,k) 
     33#   define  fse3vw_0(i,j,k)  e3vw(i,j,k) 
     34#endif 
     35#if defined key_vvl 
     36! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._1) 
     37#   define  fsdept(i,j,k)  gdept_1(i,j,k) 
     38#   define  fsdepw(i,j,k)  gdepw_1(i,j,k) 
     39#   define  fsde3w(i,j,k)  gdep3w_1(i,j,k) 
     40#   define  fse3t(i,j,k)   e3t_1(i,j,k) 
     41#   define  fse3u(i,j,k)   e3u_1(i,j,k) 
     42#   define  fse3v(i,j,k)   e3v_1(i,j,k) 
     43#   define  fse3f(i,j,k)   e3f_1(i,j,k) 
     44#   define  fse3w(i,j,k)   e3w_1(i,j,k) 
     45#   define  fse3uw(i,j,k)  e3uw_1(i,j,k) 
     46#   define  fse3vw(i,j,k)  e3vw_1(i,j,k) 
     47 
     48#   define  fsdept_b(i,j,k)  (fsdept_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))) 
     49#   define  fsdepw_b(i,j,k)  (fsdepw_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))) 
     50#   define  fsde3w_b(i,j,k)  (fsde3w_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))-sshb(i,j)) 
     51#   define  fse3t_b(i,j,k)   (fse3t_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))) 
     52#   define  fse3u_b(i,j,k)   (fse3u_0(i,j,k)*(1+sshu_b(i,j)*muu(i,j,k))) 
     53#   define  fse3v_b(i,j,k)   (fse3v_0(i,j,k)*(1+sshv_b(i,j)*muv(i,j,k))) 
     54#   define  fse3f_b(i,j,k)   (fse3f_0(i,j,k)*(1+sshf_b(i,j)*muf(i,j,k))) 
     55#   define  fse3w_b(i,j,k)   (fse3w_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))) 
     56#   define  fse3uw_b(i,j,k)  (fse3uw_0(i,j,k)*(1+sshu_b(i,j)*muu(i,j,k))) 
     57#   define  fse3vw_b(i,j,k)  (fse3vw_0(i,j,k)*(1+sshv_b(i,j)*muv(i,j,k))) 
     58 
     59#   define  fsdept_n(i,j,k)  (fsdept_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))) 
     60#   define  fsdepw_n(i,j,k)  (fsdepw_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))) 
     61#   define  fsde3w_n(i,j,k)  (fsde3w_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))-sshn(i,j)) 
     62#   define  fse3t_n(i,j,k)   (fse3t_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))) 
     63#   define  fse3u_n(i,j,k)   (fse3u_0(i,j,k)*(1+sshu_n(i,j)*muu(i,j,k))) 
     64#   define  fse3v_n(i,j,k)   (fse3v_0(i,j,k)*(1+sshv_n(i,j)*muv(i,j,k))) 
     65#   define  fse3f_n(i,j,k)   (fse3f_0(i,j,k)*(1+sshf_n(i,j)*muf(i,j,k))) 
     66#   define  fse3w_n(i,j,k)   (fse3w_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))) 
     67#   define  fse3uw_n(i,j,k)  (fse3uw_0(i,j,k)*(1+sshu_n(i,j)*muu(i,j,k))) 
     68#   define  fse3vw_n(i,j,k)  (fse3vw_0(i,j,k)*(1+sshv_n(i,j)*muv(i,j,k))) 
     69 
     70#   define  fsdept_a(i,j,k)  (fsdept_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))) 
     71#   define  fsdepw_a(i,j,k)  (fsdepw_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))) 
     72#   define  fsde3w_a(i,j,k)  (fsde3w_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))-ssha(i,j)) 
     73#   define  fse3t_a(i,j,k)   (fse3t_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))) 
     74#   define  fse3u_a(i,j,k)   (fse3u_0(i,j,k)*(1+sshu_a(i,j)*muu(i,j,k))) 
     75#   define  fse3v_a(i,j,k)   (fse3v_0(i,j,k)*(1+sshv_a(i,j)*muv(i,j,k))) 
     76#   define  fse3f_a(i,j,k)   (fse3f_0(i,j,k)*(1+sshf_a(i,j)*muf(i,j,k))) 
     77#   define  fse3w_a(i,j,k)   (fse3w_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))) 
     78#   define  fse3uw_a(i,j,k)  (fse3uw_0(i,j,k)*(1+sshu_a(i,j)*muu(i,j,k))) 
     79#   define  fse3vw_a(i,j,k)  (fse3vw_0(i,j,k)*(1+sshv_a(i,j)*muv(i,j,k))) 
     80 
     81#else 
     82! z- or s-coordinate (1D or 3D + no time dependency) use reference in all cases 
     83#   define  fsdept(i,j,k)  fsdept_0(i,j,k) 
     84#   define  fsdepw(i,j,k)  fsdepw_0(i,j,k) 
     85#   define  fsde3w(i,j,k)  fsde3w_0(i,j,k) 
     86#   define  fse3t(i,j,k)   fse3t_0(i,j,k) 
     87#   define  fse3u(i,j,k)   fse3u_0(i,j,k) 
     88#   define  fse3v(i,j,k)   fse3v_0(i,j,k) 
     89#   define  fse3f(i,j,k)   fse3f_0(i,j,k) 
     90#   define  fse3w(i,j,k)   fse3w_0(i,j,k) 
     91#   define  fse3uw(i,j,k)  fse3uw_0(i,j,k) 
     92#   define  fse3vw(i,j,k)  fse3vw_0(i,j,k) 
     93 
     94#   define  fsdept_b(i,j,k)  fsdept_0(i,j,k) 
     95#   define  fsdepw_b(i,j,k)  fsdepw_0(i,j,k) 
     96#   define  fsde3w_b(i,j,k)  fsde3w_0(i,j,k) 
     97#   define  fse3t_b(i,j,k)   fse3t_0(i,j,k) 
     98#   define  fse3u_b(i,j,k)   fse3u_0(i,j,k) 
     99#   define  fse3v_b(i,j,k)   fse3v_0(i,j,k) 
     100#   define  fse3f_b(i,j,k)   fse3f_0(i,j,k) 
     101#   define  fse3w_b(i,j,k)   fse3w_0(i,j,k) 
     102#   define  fse3uw_b(i,j,k)  fse3uw_0(i,j,k) 
     103#   define  fse3vw_b(i,j,k)  fse3vw_0(i,j,k) 
     104 
     105#   define  fsdept_n(i,j,k)  fsdept_0(i,j,k) 
     106#   define  fsdepw_n(i,j,k)  fsdepw_0(i,j,k) 
     107#   define  fsde3w_n(i,j,k)  fsde3w_0(i,j,k) 
     108#   define  fse3t_n(i,j,k)   fse3t_0(i,j,k) 
     109#   define  fse3u_n(i,j,k)   fse3u_0(i,j,k) 
     110#   define  fse3v_n(i,j,k)   fse3v_0(i,j,k) 
     111#   define  fse3f_n(i,j,k)   fse3f_0(i,j,k) 
     112#   define  fse3w_n(i,j,k)   fse3w_0(i,j,k) 
     113#   define  fse3uw_n(i,j,k)  fse3uw_0(i,j,k) 
     114#   define  fse3vw_n(i,j,k)  fse3vw_0(i,j,k) 
     115 
     116#   define  fsdept_a(i,j,k)  fsdept_0(i,j,k) 
     117#   define  fsdepw_a(i,j,k)  fsdepw_0(i,j,k) 
     118#   define  fsde3w_a(i,j,k)  fsde3w_0(i,j,k) 
     119#   define  fse3t_a(i,j,k)   fse3t_0(i,j,k) 
     120#   define  fse3u_a(i,j,k)   fse3u_0(i,j,k) 
     121#   define  fse3v_a(i,j,k)   fse3v_0(i,j,k) 
     122#   define  fse3f_a(i,j,k)   fse3f_0(i,j,k) 
     123#   define  fse3w_a(i,j,k)   fse3w_0(i,j,k) 
     124#   define  fse3uw_a(i,j,k)  fse3uw_0(i,j,k) 
     125#   define  fse3vw_a(i,j,k)  fse3vw_0(i,j,k) 
     126#endif 
     127   !!---------------------------------------------------------------------- 
     128   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    9129   !! $Id$ 
    10130   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    11    !! 
    12    !! History : 
    13    !!   9.0  !  05-10  (A. Beckmann, G. Madec) generalisation to all coord. 
    14131   !!---------------------------------------------------------------------- 
    15 #if defined key_zco 
    16 #   define  fsdept(i,j,k)  gdept_0(k) 
    17  
    18 #   define  fsdepw(i,j,k)  gdepw_0(k) 
    19 #   define  fsde3w(i,j,k)  gdepw_0(k) 
    20  
    21 #   define  fse3t(i,j,k)   e3t_0(k) 
    22 #   define  fse3u(i,j,k)   e3t_0(k) 
    23 #   define  fse3v(i,j,k)   e3t_0(k) 
    24 #   define  fse3f(i,j,k)   e3t_0(k) 
    25  
    26 #   define  fse3w(i,j,k)   e3w_0(k) 
    27 #   define  fse3uw(i,j,k)  e3w_0(k) 
    28 #   define  fse3vw(i,j,k)  e3w_0(k) 
    29 #else 
    30 #   define  fsdept(i,j,k)  gdept(i,j,k) 
    31  
    32 #   define  fsdepw(i,j,k)  gdepw(i,j,k) 
    33 #   define  fsde3w(i,j,k)  gdep3w(i,j,k) 
    34   
    35 #   define  fse3t(i,j,k)   e3t(i,j,k) 
    36 #   define  fse3u(i,j,k)   e3u(i,j,k) 
    37 #   define  fse3v(i,j,k)   e3v(i,j,k) 
    38 #   define  fse3f(i,j,k)   e3f(i,j,k) 
    39  
    40 #   define  fse3w(i,j,k)   e3w(i,j,k) 
    41 #   define  fse3uw(i,j,k)  e3uw(i,j,k) 
    42 #   define  fse3vw(i,j,k)  e3vw(i,j,k) 
    43 #endif 
    44  
    45 #if defined key_vvl 
    46 #   define  fsvdept(i,j,k)  gdept_1(i,j,k) 
    47  
    48 #   define  fsvdepw(i,j,k)  gdepw_1(i,j,k) 
    49 #   define  fsvde3w(i,j,k)  gdep3w_1(i,j,k) 
    50  
    51 #   define  fsve3t(i,j,k)   e3t_1(i,j,k) 
    52 #   define  fsve3u(i,j,k)   e3u_1(i,j,k) 
    53 #   define  fsve3v(i,j,k)   e3v_1(i,j,k) 
    54 #   define  fsve3f(i,j,k)   e3f_1(i,j,k) 
    55  
    56 #   define  fsve3w(i,j,k)   e3w_1(i,j,k) 
    57 #   define  fsve3uw(i,j,k)  e3uw_1(i,j,k) 
    58 #   define  fsve3vw(i,j,k)  e3vw_1(i,j,k) 
    59 #else 
    60 #   define  fsvdept(i,j,k)  fsdept(i,j,k) 
    61  
    62 #   define  fsvdepw(i,j,k)  fsdepw(i,j,k) 
    63 #   define  fsvde3w(i,j,k)  fsde3w(i,j,k) 
    64  
    65 #   define  fsve3t(i,j,k)   fse3t(i,j,k) 
    66 #   define  fsve3u(i,j,k)   fse3u(i,j,k) 
    67 #   define  fsve3v(i,j,k)   fse3v(i,j,k) 
    68 #   define  fsve3f(i,j,k)   fse3f(i,j,k) 
    69  
    70 #   define  fsve3w(i,j,k)   fse3w(i,j,k) 
    71 #   define  fsve3uw(i,j,k)  fse3uw(i,j,k) 
    72 #   define  fsve3vw(i,j,k)  fse3vw(i,j,k) 
    73 #endif 
    74  

  • trunk/NEMO/OPA_SRC/DYN/dynnxt.F90

    r1152 r1438  
    44   !! Ocean dynamics: time stepping 
    55   !!====================================================================== 
    6     
     6   !!====================================================================== 
     7   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code 
     8   !!                 !  1990-10  (C. Levy, G. Madec) 
     9   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions 
     10   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0 
     11   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step 
     12   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine 
     13   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module 
     14   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond. 
     15   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
     16   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.  
     17   !!            3.2  !  2009-04  (G. Madec, R.Benshila))  re-introduce the vvl option 
     18   !!---------------------------------------------------------------------- 
     19   
    720   !!---------------------------------------------------------------------- 
    821   !!   dyn_nxt      : update the horizontal velocity from the momentum trend 
    922   !!---------------------------------------------------------------------- 
    10    !! * Modules used 
    1123   USE oce             ! ocean dynamics and tracers 
    1224   USE dom_oce         ! ocean space and time domain 
    … …  
    2941   PRIVATE 
    3042 
    31    !! * Accessibility 
    32    PUBLIC dyn_nxt                ! routine called by step.F90 
     43   PUBLIC    dyn_nxt   ! routine called by step.F90 
     44 
    3345   !! * Substitutions 
    3446#  include "domzgr_substitute.h90" 
    35    !!---------------------------------------------------------------------- 
     47   !!------------------------------------------------------------------------- 
     48   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     49   !! $Id$  
     50   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
     51   !!------------------------------------------------------------------------- 
    3652 
    3753CONTAINS 
    … …  
    5975      !!             - Update un,vn arrays, the now horizontal velocity 
    6076      !! 
    61       !! History : 
    62       !!        !  87-02  (P. Andrich, D. L Hostis)  Original code 
    63       !!        !  90-10  (C. Levy, G. Madec) 
    64       !!        !  93-03  (M. Guyon)  symetrical conditions 
    65       !!        !  97-02  (G. Madec & M. Imbard)  opa, release 8.0 
    66       !!        !  97-04  (A. Weaver)  Euler forward step 
    67       !!        !  97-06  (G. Madec)  lateral boudary cond., lbc routine 
    68       !!   8.5  !  02-08  (G. Madec)  F90: Free form and module 
    69       !!        !  02-10  (C. Talandier, A-M. Treguier) Open boundary cond. 
    70       !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
    71       !!    "   !  07-07  (D. Storkey) Calls to BDY routines.  
    7277      !!---------------------------------------------------------------------- 
    73       !! * Arguments 
    7478      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    75  
    76       !! * Local declarations 
     79      !! 
    7780      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    7881      REAL(wp) ::   z2dt         ! temporary scalar 
    79       REAL(wp) ::   zsshun1, zsshvn1         ! temporary scalar 
    80       !! Variable volume 
    81       REAL(wp), DIMENSION(jpi,jpj)     ::  & ! 2D workspace 
    82          zsshub, zsshun, zsshua,           & 
    83          zsshvb, zsshvn, zsshva 
    84       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  & 
    85          zfse3ub, zfse3un, zfse3ua, &        ! 3D workspace 
    86          zfse3vb, zfse3vn, zfse3va 
    87       !!---------------------------------------------------------------------- 
    88       !!  OPA 9.0 , LOCEAN-IPSL (2005)  
    89       !! $Id$  
    90       !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    91       !!---------------------------------------------------------------------- 
     82      REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar 
     83      REAL(wp) ::   zve3a , zve3n , zve3b    !    -         - 
     84      REAL(wp) ::   ze3u_b, ze3u_n, ze3u_a   !    -         - 
     85      REAL(wp) ::   ze3v_b, ze3v_n, ze3v_a   !    -         -  
     86      REAL(wp) ::   zuf   , zvf              !    -         -  
     87     !!---------------------------------------------------------------------- 
    9288 
    9389      IF( kt == nit000 ) THEN 
    … …  
    10298 
    10399      !! Explicit physics with thickness weighted updates 
    104       IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN 
    105  
    106          ! Sea surface elevation time stepping 
    107          ! ----------------------------------- 
    108          ! 
    109          DO jj = 1, jpjm1 
    110             DO ji = 1,jpim1 
    111  
    112                ! Sea Surface Height at u-point before 
    113                zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )     & 
    114                   &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * sshbb(ji  ,jj  )   & 
    115                   &            + e1t(ji+1,jj  ) * e2t(ji+1,jj  ) * sshbb(ji+1,jj  ) ) 
    116  
    117                ! Sea Surface Height at v-point before 
    118                zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )     & 
    119                   &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * sshbb(ji  ,jj  )   & 
    120                   &            + e1t(ji  ,jj+1) * e2t(ji  ,jj+1) * sshbb(ji  ,jj+1) ) 
    121  
    122                ! Sea Surface Height at u-point after 
    123                zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )     & 
    124                   &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * ssha(ji  ,jj  )    & 
    125                   &            + e1t(ji+1,jj  ) * e2t(ji+1,jj  ) * ssha(ji+1,jj  ) ) 
    126  
    127                ! Sea Surface Height at v-point after 
    128                zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )     & 
    129                   &          * ( e1t(ji  ,jj  ) * e2t(ji  ,jj  ) * ssha(ji  ,jj  )    & 
    130                   &            + e1t(ji  ,jj+1) * e2t(ji  ,jj+1) * ssha(ji  ,jj+1) ) 
    131  
    132             END DO 
    133          END DO 
    134  
    135          ! Boundaries conditions 
    136          CALL lbc_lnk( zsshua, 'U', 1. )   ;    CALL lbc_lnk( zsshva, 'V', 1. ) 
    137          CALL lbc_lnk( zsshub, 'U', 1. )   ;    CALL lbc_lnk( zsshvb, 'V', 1. ) 
    138  
    139          ! Scale factors at before and after time step 
    140          ! ------------------------------------------- 
    141          CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ;    CALL dom_vvl_sf( zsshua, 'U', zfse3ua )  
    142          CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ;    CALL dom_vvl_sf( zsshva, 'V', zfse3va )  
    143  
    144          ! Asselin filtered scale factor at now time step 
    145          ! ---------------------------------------------- 
    146          IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN 
    147             CALL dom_vvl_sf_ini( 'U', zfse3un ) ;   CALL dom_vvl_sf_ini( 'V', zfse3vn )  
    148          ELSE 
    149             zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:) 
    150             zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:) 
    151             CALL dom_vvl_sf( zsshun, 'U', zfse3un ) ;   CALL dom_vvl_sf( zsshvn, 'V', zfse3vn )  
    152          ENDIF 
    153  
    154          ! Thickness weighting 
    155          ! ------------------- 
     100 
     101      ! Lateral boundary conditions on ( ua, va ) 
     102      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )  
     103 
     104      ! Next velocity 
     105      ! ------------- 
     106#if defined key_dynspg_flt 
     107      ! Leap-frog time stepping already done in dynspg_flt.F routine 
     108#else 
     109      IF( lk_vvl ) THEN                                  ! Varying levels 
    156110         DO jk = 1, jpkm1 
    157             DO jj = 1, jpj 
    158                DO ji = 1, jpi 
    159                   ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk) 
    160                   va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk) 
    161  
    162                   un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk) 
    163                   vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk) 
    164  
    165                   ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk) 
    166                   vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk) 
    167                END DO 
    168             END DO 
    169          END DO 
    170  
    171       ENDIF 
    172  
    173       ! Lateral boundary conditions on ( ua, va ) 
    174       CALL lbc_lnk( ua, 'U', -1. ) 
    175       CALL lbc_lnk( va, 'V', -1. ) 
    176  
    177       !                                                ! =============== 
    178       DO jk = 1, jpkm1                                 ! Horizontal slab 
    179          !                                             ! =============== 
    180          ! Next velocity 
    181          ! ------------- 
    182 #if defined key_dynspg_flt 
    183          ! Leap-frog time stepping already done in dynspg.F routine 
    184 #else 
    185          DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
    186             DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    187                ! Leap-frog time stepping 
    188                ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 
    189                va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 
    190             END DO 
    191          END DO 
    192  
    193          IF( lk_vvl ) THEN 
    194             ! Unweight velocities prior to updating open boundaries. 
    195  
    196             DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
    197                DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    198                   ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk) 
    199                   va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk) 
    200  
    201                   un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk) 
    202                   vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk) 
    203  
    204                   ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk) 
    205                   vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk) 
    206                END DO 
    207             END DO 
    208  
    209          ENDIF 
     111            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)     & 
     112               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) )   & 
     113               &         / fse3u_a(:,:,jk) * umask(:,:,jk) 
     114            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)     & 
     115               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) )   & 
     116               &         / fse3v_a(:,:,jk) * vmask(:,:,jk) 
     117         END DO 
     118      ELSE 
     119         DO jk = 1, jpkm1 
     120                  ! Leap-frog time stepping 
     121            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     122            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     123         END DO 
     124      ENDIF 
    210125 
    211126# if defined key_obc 
    212          !                                             ! =============== 
    213       END DO                                           !   End of slab 
    214       !                                                ! =============== 
    215127      ! Update (ua,va) along open boundaries (only in the rigid-lid case) 
    216128      CALL obc_dyn( kt ) 
    217129 
    218130      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 
    219          !Flather boundary condition : 
     131         ! Flather boundary condition : 
    220132         !        - Update sea surface height on each open boundary 
    221133         !                 sshn (= after ssh) for explicit case 
    … …  
    223135         !        - Correct the barotropic velocities 
    224136         CALL obc_dyn_bt( kt ) 
    225  
    226          !Boundary conditions on sshn ( after ssh) 
     137         ! 
     138         ! Boundary conditions on sshn ( after ssh) 
    227139         CALL lbc_lnk( sshn, 'T', 1. ) 
    228  
    229          IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    230             CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask) 
    231          ENDIF 
    232  
    233          IF ( ln_vol_cst ) CALL obc_vol( kt ) 
    234  
    235       ENDIF 
    236  
    237       !                                                ! =============== 
    238       DO jk = 1, jpkm1                                 ! Horizontal slab 
    239          !                                             ! =============== 
     140         ! 
     141         IF( ln_vol_cst )   CALL obc_vol( kt ) 
     142         ! 
     143         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask ) 
     144      ENDIF 
     145 
    240146# elif defined key_bdy  
    241          !                                             ! =============== 
    242       END DO                                           !   End of slab 
    243       !                                                ! =============== 
    244147      ! Update (ua,va) along open boundaries (for exp or ts options). 
    245       IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN 
    246    
     148      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 
     149         ! 
    247150         CALL bdy_dyn_frs( kt ) 
    248  
    249          IF ( ln_bdy_fla ) THEN 
    250  
    251             ua_e(:,:)=0.0 
    252             va_e(:,:)=0.0 
    253  
     151         ! 
     152         IF( ln_bdy_fla ) THEN 
     153            ua_e(:,:) = 0.e0 
     154            va_e(:,:) = 0.e0 
    254155            ! Set these variables for use in bdy_dyn_fla 
    255156            hu_e(:,:) = hu(:,:) 
    256157            hv_e(:,:) = hv(:,:) 
    257  
    258             DO jk = 1, jpkm1 
    259                !! Vertically integrated momentum trends 
     158            DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    260159               ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    261160               va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    262161            END DO 
    263  
    264162            DO jk = 1 , jpkm1 
    265163               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:) 
    266164               va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:) 
    267165            END DO 
    268  
    269166            CALL bdy_dta_bt( kt+1, 0) 
    270167            CALL bdy_dyn_fla 
    271  
    272168         ENDIF 
    273  
     169         ! 
    274170         DO jk = 1 , jpkm1 
    275171            ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:) 
    276172            va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:) 
    277173         END DO 
    278  
    279       ENDIF 
    280  
    281       !                                                ! =============== 
    282       DO jk = 1, jpkm1                                 ! Horizontal slab 
    283          !                                             ! =============== 
     174      ENDIF 
    284175# endif 
     176 
    285177# if defined key_agrif 
    286          !                                             ! =============== 
    287       END DO                                           !   End of slab 
    288       !                                                ! =============== 
    289       ! Update (ua,va) along open boundaries (only in the rigid-lid case) 
    290178      CALL Agrif_dyn( kt ) 
    291       !                                                ! =============== 
    292       DO jk = 1, jpkm1                                 ! Horizontal slab 
    293          !                                             ! =============== 
    294179# endif 
    295180#endif 
    296181 
    297          ! Time filter and swap of dynamics arrays 
    298          ! ------------------------------------------ 
    299          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    300             IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels 
    301                ! caution: don't use (:,:) for this loop 
    302                ! it causes optimization problems on NEC in auto-tasking 
     182      ! Time filter and swap of dynamics arrays 
     183      ! ------------------------------------------ 
     184      IF( neuler == 0 .AND. kt == nit000 ) THEN 
     185         DO jk = 1, jpkm1                                  
     186            un(:,:,jk) = ua(:,:,jk) 
     187            vn(:,:,jk) = va(:,:,jk) 
     188         END DO 
     189      ELSE 
     190         IF( lk_vvl ) THEN                                ! Varying levels 
     191            DO jk = 1, jpkm1                              ! filter applied on thickness weighted velocities 
    303192               DO jj = 1, jpj 
    304193                  DO ji = 1, jpi 
    305                      zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk) 
    306                      zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk) 
    307                      ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk) 
    308                      vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk) 
    309                      zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 
    310                      zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 
    311                      un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk) 
    312                      vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk) 
     194                     ze3u_a = fse3u_a(ji,jj,jk) 
     195                     ze3v_a = fse3v_a(ji,jj,jk) 
     196                     ze3u_n = fse3u_n(ji,jj,jk) 
     197                     ze3v_n = fse3v_n(ji,jj,jk) 
     198                     ze3u_b = fse3u_b(ji,jj,jk) 
     199                     ze3v_b = fse3v_b(ji,jj,jk) 
     200                     ! 
     201                     zue3a = ua(ji,jj,jk) * ze3u_a 
     202                     zve3a = va(ji,jj,jk) * ze3v_a 
     203                     zue3n = un(ji,jj,jk) * ze3u_n 
     204                     zve3n = vn(ji,jj,jk) * ze3v_n 
     205                     zue3b = ub(ji,jj,jk) * ze3u_b 
     206                     zve3b = vb(ji,jj,jk) * ze3v_b 
     207                     ! 
     208                     ub(ji,jj,jk) = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   & 
     209                        &         / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 
     210                     vb(ji,jj,jk) = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   & 
     211                        &         / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 
     212                     un(ji,jj,jk) = ua(ji,jj,jk) * umask(ji,jj,jk) 
     213                     vn(ji,jj,jk) = va(ji,jj,jk) * vmask(ji,jj,jk) 
    313214                  END DO 
    314215               END DO 
    315             ELSE                                               ! Fixed levels 
     216            END DO 
     217         ELSE                                             ! Fixed levels 
     218            DO jk = 1, jpkm1                              ! filter applied on velocities 
    316219               DO jj = 1, jpj 
    317220                  DO ji = 1, jpi 
    318                      ! Euler (forward) time stepping 
    319                      ub(ji,jj,jk) = un(ji,jj,jk) 
    320                      vb(ji,jj,jk) = vn(ji,jj,jk) 
     221                     zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 
     222                     zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 
     223                     ub(ji,jj,jk) = zuf 
     224                     vb(ji,jj,jk) = zvf 
    321225                     un(ji,jj,jk) = ua(ji,jj,jk) 
    322226                     vn(ji,jj,jk) = va(ji,jj,jk) 
    323227                  END DO 
    324228               END DO 
    325             ENDIF 
    326          ELSE 
    327             IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels 
    328                ! caution: don't use (:,:) for this loop 
    329                ! it causes optimization problems on NEC in auto-tasking 
    330                DO jj = 1, jpj 
    331                   DO ji = 1, jpi 
    332                      zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk) 
    333                      zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk) 
    334                      ub(ji,jj,jk) = ( atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) & 
    335                        &            + atfp1 * un(ji,jj,jk) ) * zsshun1 
    336                      vb(ji,jj,jk) = ( atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) & 
    337                        &            + atfp1 * vn(ji,jj,jk) ) * zsshvn1 
    338                      zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 
    339                      zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 
    340                      un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 
    341                      vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 
    342                   END DO 
    343                END DO 
    344             ELSE                                               ! Fixed levels 
    345                DO jj = 1, jpj 
    346                   DO ji = 1, jpi 
    347                      ! Leap-frog time stepping 
    348                      ub(ji,jj,jk) = atfp  * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 
    349                      vb(ji,jj,jk) = atfp  * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 
    350                      un(ji,jj,jk) = ua(ji,jj,jk) 
    351                      vn(ji,jj,jk) = va(ji,jj,jk) 
    352                   END DO 
    353                END DO 
    354             ENDIF 
     229            END DO 
    355230         ENDIF 
    356          !                                             ! =============== 
    357       END DO                                           !   End of slab 
    358       !                                                ! =============== 
    359  
    360       IF(ln_ctl)   THEN 
    361          CALL prt_ctl(tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask, & 
    362             &         tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask) 
    363231      ENDIF 
    364232 
    … …  
    367235#endif       
    368236 
     237      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
     238         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
     239      !  
    369240   END SUBROUTINE dyn_nxt 
    370241 

  • trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90

    r1146 r1438  
    44   !! Ocean dynamics:  surface pressure gradient trend 
    55   !!====================================================================== 
     6   !! History :  9.0  !  2005-11  (V. Garnier, G. Madec, L. Bessieres) Original code 
     7   !!---------------------------------------------------------------------- 
    68#if defined key_dynspg_exp   ||   defined key_esopa 
    79   !!---------------------------------------------------------------------- 
    … …  
    1113   !!                      pressure gradient in the free surface constant   
    1214   !!                      volume case with vector optimization 
    13    !!   exp_rst      : read/write the explicit restart fields in the ocean restart file 
    1415   !!---------------------------------------------------------------------- 
    15    !! * Modules used 
    1616   USE oce             ! ocean dynamics and tracers  
    1717   USE dom_oce         ! ocean space and time domain  
    … …  
    3131   PRIVATE 
    3232 
    33    !! * Accessibility 
    34    PUBLIC dyn_spg_exp  ! routine called by step.F90 
    35    PUBLIC exp_rst      ! routine called by istate.F90 
     33   PUBLIC   dyn_spg_exp   ! routine called by step.F90 
    3634 
    3735   !! * Substitutions 
    … …  
    3937#  include "vectopt_loop_substitute.h90" 
    4038   !!---------------------------------------------------------------------- 
    41    !!  OPA 9.0 , LOCEAN-IPSL (2005)  
     39   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    4240   !! $Id$ 
    43    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     41   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
    4442   !!---------------------------------------------------------------------- 
    45  
    4643 
    4744CONTAINS 
    … …  
    6259      !!      -1- Compute the now surface pressure gradient 
    6360      !!      -2- Add it to the general trend 
    64       !!      -3- Compute the horizontal divergence of velocities 
    65       !!      - the now divergence is given by : 
    66       !!         zhdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    67       !!      - integrate the horizontal divergence from the bottom  
    68       !!         to the surface 
    69       !!      - apply lateral boundary conditions on zhdivn 
    70       !!      -4- Estimate the after sea surface elevation from the kinematic 
    71       !!         surface boundary condition: 
    72       !!         zssha = sshb - 2 rdt ( zhdiv + emp ) 
    73       !!      - Time filter applied on now sea surface elevation to avoid 
    74       !!         the divergence of two consecutive time-steps and swap of free 
    75       !!         surface arrays to start the next time step: 
    76       !!         sshb = sshn + atfp * [ sshb + zssha - 2 sshn ] 
    77       !!         sshn = zssha 
    78       !!      - apply lateral boundary conditions on sshn 
    7961      !! 
    8062      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 
     63      !!--------------------------------------------------------------------- 
     64      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index 
    8165      !! 
    82       !! References : 
    83       !! 
    84       !! History : 
    85       !!   9.0  !  05-11  (V. Garnier, G. Madec, L. Bessieres) Original code 
    86       !! 
    87       !!--------------------------------------------------------------------- 
    88       !! * Arguments 
    89       INTEGER, INTENT( in )  ::   kt         ! ocean time-step index 
    90  
    91       !! * Local declarations 
    9266      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
    93       REAL(wp) ::   z2dt, zraur, zssha       ! temporary scalars  
    94       REAL(wp), DIMENSION(jpi,jpj)    ::  &  ! temporary arrays 
    95          &         zhdiv 
    9667      !!---------------------------------------------------------------------- 
    9768 
    … …  
    10475         spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction) 
    10576         spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction) 
    106  
    107          CALL exp_rst( nit000, 'READ' )       ! read or initialize the following fields: 
    108          !                                    ! sshb, sshn 
    109  
    11077      ENDIF 
    11178 
    112       IF( .NOT. lk_vvl ) THEN 
     79      ! read or estimate sea surface height and vertically integrated velocities 
     80      IF( lk_obc )   CALL obc_dta_bt( kt, 0 ) 
    11381 
    114          ! 0. Initialization 
    115          ! ----------------- 
    116          ! read or estimate sea surface height and vertically integrated velocities 
    117          IF( lk_obc )   CALL obc_dta_bt( kt, 0 ) 
    118          z2dt = 2. * rdt                                       ! time step: leap-frog 
    119          IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest 
    120          zraur = 1.e0 / rauw 
     82      ! Surface pressure gradient (now) 
     83      DO jj = 2, jpjm1 
     84         DO ji = fs_2, fs_jpim1   ! vector opt. 
     85            spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 
     86            spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 
     87         END DO 
     88      END DO 
    12189 
    122          ! 1. Surface pressure gradient (now) 
    123          ! ---------------------------- 
     90      ! Add the surface pressure trend to the general trend 
     91      DO jk = 1, jpkm1 
    12492         DO jj = 2, jpjm1 
    12593            DO ji = fs_2, fs_jpim1   ! vector opt. 
    126                spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 
    127                spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 
     94               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 
     95               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 
    12896            END DO 
    12997         END DO 
     98      END DO 
    13099 
    131          ! 2. Add the surface pressure trend to the general trend 
    132          ! ------------------------------------------------------ 
    133          DO jk = 1, jpkm1 
    134             DO jj = 2, jpjm1 
    135                DO ji = fs_2, fs_jpim1   ! vector opt. 
    136                   ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 
    137                   va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 
    138                END DO 
    139             END DO 
    140          END DO 
    141       
    142          ! 3. Vertical integration of horizontal divergence of velocities 
    143          ! -------------------------------- 
    144          zhdiv(:,:) = 0.e0 
    145          DO jk = jpkm1, 1, -1 
    146             DO jj = 2, jpjm1 
    147                DO ji = fs_2, fs_jpim1   ! vector opt. 
    148                   zhdiv(ji,jj) = zhdiv(ji,jj) + (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)      & 
    149                      &                           - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)      & 
    150                      &                           + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)      & 
    151                      &                           - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )   & 
    152                      &                        / ( e1t(ji,jj) * e2t(ji,jj) ) 
    153                END DO 
    154             END DO 
    155          END DO 
    156  
    157 #if defined key_obc 
    158          ! open boundaries (div must be zero behind the open boundary) 
    159          !  mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    160          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east 
    161          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west 
    162          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north 
    163          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south 
    164 #endif 
    165  
    166          ! 4. Sea surface elevation time stepping 
    167          ! -------------------------------------- 
    168          ! Euler (forward) time stepping, no time filter 
    169          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    170             DO jj = 1, jpj 
    171                DO ji = 1, jpi 
    172                   ! after free surface elevation 
    173                   zssha = sshb(ji,jj) - rdt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1) 
    174                   ! swap of arrays 
    175                   sshb(ji,jj) = sshn(ji,jj) 
    176                   sshn(ji,jj) = zssha 
    177                END DO 
    178             END DO 
    179          ELSE 
    180             ! Leap-frog time stepping and time filter 
    181             DO jj = 1, jpj 
    182                DO ji = 1, jpi 
    183                   ! after free surface elevation 
    184                   zssha = sshb(ji,jj) - z2dt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1) 
    185                   ! time filter and array swap 
    186                   sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj) 
    187                   sshn(ji,jj) = zssha 
    188                END DO 
    189             END DO 
    190          ENDIF 
    191  
    192       ELSE !! Variable volume, ssh time-stepping already done 
    193  
    194          ! read or estimate sea surface height and vertically integrated velocities 
    195          IF( lk_obc )   CALL obc_dta_bt( kt, 0 ) 
    196  
    197          ! Sea surface elevation swap 
    198          ! ----------------------------- 
    199          ! 
    200          sshbb(:,:) = sshb(:,:) ! Needed for the dynamics time-stepping 
    201          ! 
    202          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    203             ! No time filter 
    204             sshb(:,:) = sshn(:,:) 
    205             sshn(:,:) = ssha(:,:) 
    206          ELSE 
    207             ! Time filter 
    208             sshb(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:) 
    209             sshn(:,:) = ssha(:,:) 
    210          ENDIF 
    211  
    212       ENDIF 
    213  
    214       ! Boundary conditions on sshn 
    215       IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. ) 
    216  
    217       ! write filtered free surface arrays in restart file 
    218       ! -------------------------------------------------- 
    219       IF( lrst_oce )   CALL exp_rst( kt, 'WRITE' ) 
    220   
    221       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    222          CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask) 
    223       ENDIF 
    224        
    225100   END SUBROUTINE dyn_spg_exp 
    226101 
    227    SUBROUTINE exp_rst( kt, cdrw ) 
    228       !!--------------------------------------------------------------------- 
    229       !!                   ***  ROUTINE exp_rst  *** 
    230       !! 
    231       !! ** Purpose : Read or write explicit arrays in restart file 
    232       !!---------------------------------------------------------------------- 
    233       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    234       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    235       ! 
    236       !!---------------------------------------------------------------------- 
    237       ! 
    238       IF( TRIM(cdrw) == 'READ' ) THEN 
    239          IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
    240             CALL iom_get( numror, jpdom_autoglo, 'sshb'  , sshb(:,:)   ) 
    241             CALL iom_get( numror, jpdom_autoglo, 'sshn'  , sshn(:,:)   ) 
    242             IF( neuler == 0 ) sshb(:,:) = sshn(:,:) 
    243          ELSE 
    244             IF( nn_rstssh == 1 ) THEN   
    245                sshb(:,:) = 0.e0 
    246                sshn(:,:) = 0.e0 
    247             ENDIF 
    248          ENDIF 
    249       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    250          CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) ) 
    251          CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) ) 
    252       ENDIF 
    253       ! 
    254    END SUBROUTINE exp_rst 
    255102#else 
    256103   !!---------------------------------------------------------------------- 
    … …  
    261108      WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt 
    262109   END SUBROUTINE dyn_spg_exp 
    263    SUBROUTINE exp_rst( kt, cdrw )     ! Empty routine 
    264       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    265       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    266       WRITE(*,*) 'exp_rst: You should not have seen this print! error?', kt, cdrw 
    267    END SUBROUTINE exp_rst 
    268110#endif 
    269     
     111 
    270112   !!====================================================================== 
    271113END MODULE dynspg_exp 

  • trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r1200 r1438  
    44   !! Ocean dynamics:  surface pressure gradient trend 
    55   !!====================================================================== 
    6    !! History    8.0  !  98-05  (G. Roullet)  free surface 
    7    !!                 !  98-10  (G. Madec, M. Imbard)  release 8.2 
    8    !!            8.5  !  02-08  (G. Madec)  F90: Free form and module 
    9    !!            " "  !  02-11  (C. Talandier, A-M Treguier) Open boundaries 
    10    !!            9.0  !  04-08  (C. Talandier) New trends organization 
    11    !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization 
    12    !!            " "  !  06-07  (S. Masson)  distributed restart using iom 
    13    !!            " "  !  05-01  (J.Chanut, A.Sellar) Calls to BDY routines.  
     6   !! History    OPA  !  1998-05  (G. Roullet)  free surface 
     7   !!                 !  1998-10  (G. Madec, M. Imbard)  release 8.2 
     8   !!   NEMO     O.1  !  2002-08  (G. Madec)  F90: Free form and module 
     9   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries 
     10   !!            1.0  !  2004-08  (C. Talandier) New trends organization 
     11   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization 
     12   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom 
     13   !!             -   !  2006-08  (J.Chanut, A.Sellar) Calls to BDY routines.  
     14   !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 
    1415   !!---------------------------------------------------------------------- 
    1516#if defined key_dynspg_flt   ||   defined key_esopa   
    1617   !!---------------------------------------------------------------------- 
    1718   !!   'key_dynspg_flt'                              filtered free surface 
    18    !!---------------------------------------------------------------------- 
    1919   !!---------------------------------------------------------------------- 
    2020   !!   dyn_spg_flt  : update the momentum trend with the surface pressure 
    … …  
    5353   PRIVATE 
    5454 
    55    PUBLIC dyn_spg_flt  ! routine called by step.F90 
    56    PUBLIC flt_rst      ! routine called by istate.F90 
     55   PUBLIC   dyn_spg_flt  ! routine called by step.F90 
     56   PUBLIC   flt_rst      ! routine called by istate.F90 
    5757 
    5858   !! * Substitutions 
    … …  
    6060#  include "vectopt_loop_substitute.h90" 
    6161   !!---------------------------------------------------------------------- 
    62    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
     62   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    6363   !! $Id$ 
    6464   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
    … …  
    8080      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + rnu btda ) 
    8181      !!      where sshn is the free surface elevation and btda is the after 
    82       !!      of the free surface elevation 
    83       !!       -1- compute the after sea surface elevation from the kinematic 
    84       !!      surface boundary condition: 
    85       !!              zssha = sshb + 2 rdt ( wn - emp ) 
    86       !!           Time filter applied on now sea surface elevation to avoid  
    87       !!      the divergence of two consecutive time-steps and swap of free 
    88       !!      surface arrays to start the next time step: 
    89       !!              sshb = sshn + atfp * [ sshb + zssha - 2 sshn ] 
    90       !!              sshn = zssha 
    91       !!       -2- evaluate the surface presure trend (including the addi- 
     82      !!      time derivative of the free surface elevation 
     83      !!       -1- evaluate the surface presure trend (including the addi- 
    9284      !!      tional force) in three steps: 
    9385      !!        a- compute the right hand side of the elliptic equation: 
    … …  
    10597      !!         iterative solver 
    10698      !!        c- apply the solver by a call to sol... routine 
    107       !!       -3- compute and add the free surface pressure gradient inclu- 
     99      !!       -2- compute and add the free surface pressure gradient inclu- 
    108100      !!      ding the additional force used to stabilize the equation. 
    109101      !! 
    … …  
    112104      !! References : Roullet and Madec 1999, JGR. 
    113105      !!--------------------------------------------------------------------- 
    114       !! * Modules used 
    115       USE oce    , ONLY :   zub   => ta,             &  ! ta used as workspace 
    116                             zvb   => sa                 ! sa   "          " 
    117  
    118       INTEGER, INTENT( in )  ::   kt         ! ocean time-step index 
    119       INTEGER, INTENT( out ) ::   kindic     ! solver convergence flag (<0 if not converge) 
     106      USE oce, ONLY :   zub   => ta   ! ta used as workspace 
     107      USE oce, ONLY :   zvb   => sa   ! ta used as workspace 
     108      !! 
     109      INTEGER, INTENT( in )  ::   kt       ! ocean time-step index 
     110      INTEGER, INTENT( out ) ::   kindic   ! solver convergence flag (<0 if not converge) 
    120111      !!                                    
    121112      INTEGER  ::   ji, jj, jk                          ! dummy loop indices 
    122       REAL(wp) ::   z2dt, z2dtg, zraur, znugdt,      &  ! temporary scalars 
    123          &          znurau, zgcb, zbtd,              &  !   "          " 
    124          &          ztdgu, ztdgv                        !   "          " 
    125       !! Variable volume 
    126       REAL(wp), DIMENSION(jpi,jpj)     ::  & 
    127          zsshub, zsshua, zsshvb, zsshva, zssha          ! 2D workspace 
    128       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &            ! 3D workspace 
    129          zfse3ub, zfse3ua, zfse3vb, zfse3va             ! 3D workspace 
     113      REAL(wp) ::   z2dt, z2dtg, zraur, znugdt          ! temporary scalars 
     114      REAL(wp) ::   znurau, zgcb, zbtd                  !   "          " 
     115      REAL(wp) ::   ztdgu, ztdgv                        !   "          " 
    130116      !!---------------------------------------------------------------------- 
    131117      ! 
    … …  
    143129         ! when using agrif, sshn, gcx have to be read in istate 
    144130         IF (.NOT. lk_agrif) CALL flt_rst( nit000, 'READ' )       ! read or initialize the following fields: 
    145          !                                                        ! gcx, gcxb, sshb, sshn 
     131         !                                                        ! gcx, gcxb 
    146132      ENDIF 
    147133 
    … …  
    158144      IF( lk_vvl ) THEN          ! variable volume  
    159145 
    160          DO jj = 1, jpjm1 
    161             DO ji = 1,jpim1 
    162  
    163                ! Sea Surface Height at u-point before 
    164                zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  &  
    165                   &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)       & 
    166                   &            + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
    167  
    168                ! Sea Surface Height at v-point before 
    169                zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  & 
    170                   &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )       & 
    171                   &            + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    172  
    173                ! Sea Surface Height at u-point after 
    174                zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  & 
    175                   &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)       & 
    176                   &            + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
    177  
    178                ! Sea Surface Height at v-point after 
    179                zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  & 
    180                   &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )       & 
    181                   &            + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
    182  
    183             END DO 
    184          END DO 
    185  
    186          ! Boundaries conditions 
    187          CALL lbc_lnk( zsshub, 'U', 1. )    ;    CALL lbc_lnk( zsshvb, 'V', 1. ) 
    188          CALL lbc_lnk( zsshua, 'U', 1. )    ;    CALL lbc_lnk( zsshva, 'V', 1. ) 
    189  
    190          ! Scale factors at before and after time step 
    191          ! ------------------------------------------- 
    192          CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ;    CALL dom_vvl_sf( zsshua, 'U', zfse3ua ) 
    193          CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ;    CALL dom_vvl_sf( zsshva, 'V', zfse3va ) 
    194  
    195  
    196          ! Thickness weighting 
    197          ! ------------------- 
    198          DO jk = 1, jpkm1 
    199             DO jj = 1, jpj 
    200                DO ji = 1, jpi 
    201                   ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk) 
    202                   va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk) 
    203  
    204                   zub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk) 
    205                   zvb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk) 
    206                END DO 
    207             END DO 
    208          END DO 
    209  
    210146         ! Evaluate the masked next velocity (effect of the additional force not included) 
     147         ! -------------------   (thickness weighted velocity, surface pressure gradient already included in dyn_hpg) 
    211148         DO jk = 1, jpkm1 
    212149            DO jj = 2, jpjm1 
    213150               DO ji = fs_2, fs_jpim1   ! vector opt. 
    214                   ua(ji,jj,jk) = ( zub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) /zfse3ua(ji,jj,jk) * umask(ji,jj,jk) 
    215                   va(ji,jj,jk) = ( zvb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) /zfse3va(ji,jj,jk) * vmask(ji,jj,jk) 
     151                  ua(ji,jj,jk) = (        ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      & 
     152                     &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   & 
     153                     &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk) 
     154                  va(ji,jj,jk) = (        vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      & 
     155                     &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   & 
     156                     &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk) 
    216157               END DO 
    217158            END DO 
    … …  
    227168            END DO  
    228169         END DO  
    229  
    230          ! Add the surface pressure trend to the general trend 
     170         ! 
     171         ! add the surface pressure trend to the general trend and 
     172         ! evaluate the masked next velocity (effect of the additional force not included) 
    231173         DO jk = 1, jpkm1 
    232174            DO jj = 2, jpjm1 
    233175               DO ji = fs_2, fs_jpim1   ! vector opt. 
    234                   ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 
    235                   va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 
     176                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk) 
     177                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk) 
    236178               END DO 
    237179            END DO 
    238180         END DO 
    239  
    240          ! Evaluate the masked next velocity (effect of the additional force not included) 
    241          DO jk = 1, jpkm1 
    242             DO jj = 2, jpjm1 
    243                DO ji = fs_2, fs_jpim1   ! vector opt. 
    244                   ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 
    245                   va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 
    246                END DO 
    247             END DO 
    248          END DO 
    249  
     181         ! 
    250182      ENDIF 
    251183 
    … …  
    311243      CALL lbc_lnk( spgv, 'V', -1. ) 
    312244 
    313       IF( lk_vvl ) CALL sol_mat( kt ) 
     245      IF( lk_vvl ) CALL sol_mat( kt )      ! build the matrix at kt (vvl case only) 
    314246 
    315247      ! Right hand side of the elliptic equation and first guess 
    … …  
    347279      ! Relative precision (computation on one processor) 
    348280      ! ------------------ 
    349       rnorme =0. 
     281      rnorme =0.e0 
    350282      rnorme = SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) ) 
    351283      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain 
    … …  
    413345      ENDIF 
    414346#endif       
    415       ! 7.  Add the trends multiplied by z2dt to the after velocity 
    416       ! ----------------------------------------------------------- 
     347      ! Add the trends multiplied by z2dt to the after velocity 
     348      ! ------------------------------------------------------- 
    417349      !     ( c a u t i o n : (ua,va) here are the after velocity not the 
    418350      !                       trend, the leap-frog time stepping will not 
    … …  
    421353         DO jj = 2, jpjm1 
    422354            DO ji = fs_2, fs_jpim1   ! vector opt. 
    423                ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk) 
    424                va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk) 
     355               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk) 
     356               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk) 
    425357            END DO 
    426358         END DO 
    427359      END DO 
    428  
    429       ! Sea surface elevation time stepping 
    430       ! ----------------------------------- 
    431       IF( lk_vvl ) THEN   ! after free surface elevation 
    432          zssha(:,:) = ssha(:,:) 
    433       ELSE 
    434          zssha(:,:) = sshb(:,:) + z2dt * ( wn(:,:,1) - emp(:,:) * zraur ) * tmask(:,:,1) 
    435       ENDIF 
    436 #if defined key_obc 
    437 # if defined key_agrif 
    438       IF ( Agrif_Root() ) THEN 
    439 # endif 
    440          zssha(:,:)=zssha(:,:)*obctmsk(:,:) 
    441          CALL lbc_lnk(zssha,'T',1.)  ! absolutly compulsory !! (jmm) 
    442 # if defined key_agrif 
    443       ENDIF 
    444 # endif 
    445 #endif 
    446  
    447       IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler (forward) time stepping, no time filter 
    448          ! swap of arrays 
    449          sshb(:,:) = sshn (:,:) 
    450          sshn(:,:) = zssha(:,:) 
    451       ELSE                                           ! Leap-frog time stepping and time filter 
    452          ! time filter and array swap 
    453          sshb(:,:) = atfp * ( sshb(:,:) + zssha(:,:) ) + atfp1 * sshn(:,:) 
    454          sshn(:,:) = zssha(:,:) 
    455       ENDIF 
    456360 
    457361      ! write filtered free surface arrays in restart file 
    458362      ! -------------------------------------------------- 
    459363      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 
    460  
    461       ! print sum trends (used for debugging) 
    462       IF(ln_ctl)     CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask ) 
    463364      ! 
    464365   END SUBROUTINE dyn_spg_flt 
    … …  
    481382           CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) ) 
    482383           CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) ) 
    483            CALL iom_get( numror, jpdom_autoglo, 'sshb', sshb(:,:)         ) 
    484            CALL iom_get( numror, jpdom_autoglo, 'sshn', sshn(:,:)         ) 
    485            IF( neuler == 0 ) THEN 
    486               sshb(:,:) = sshn(:,:) 
    487               gcxb(:,:) = gcx (:,:) 
    488            ENDIF 
     384           IF( neuler == 0 )   gcxb(:,:) = gcx (:,:) 
    489385        ELSE 
    490386           gcx (:,:) = 0.e0 
    491387           gcxb(:,:) = 0.e0 
    492            IF( nn_rstssh == 1 ) THEN   
    493               sshb(:,:) = 0.e0 
    494               sshn(:,:) = 0.e0 
    495            ENDIF 
    496388        ENDIF 
    497389     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    498390! Caution : extra-hallow 
    499391! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 
    500         CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) ) 
     392        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) ) 
    501393        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 
    502         CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:)         ) 
    503         CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:)         ) 
    504394     ENDIF 
    505395     ! 

  • trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90

    r1152 r1438  
    4343   !! ------------------------ 
    4444      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop 
    45          sshn_b, sshb_b,               &  ! sea surface heigth (now, before) 
    46          un_b  , vn_b                     ! vertically integrated horizontal velocities (now) 
     45         un_e  , vn_e,                        & ! vertically integrated horizontal velocities (now) 
     46         ua_e  , va_e                           ! vertically integrated horizontal velocities (after) 
    4747      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop 
    48          sshb_e, sshn_e, ssha_e,       &  ! sea surface heigth (before,now,after) 
    49          ua_e  , va_e,                 &  ! vertically integrated horizontal velocities (after) 
    50          hu_e  , hv_e                     ! depth arrays for the barotropic solution 
     48         sshn_e, ssha_e,                      & ! sea surface heigth (now, after) 
     49         hu_e  , hv_e                           ! depth arrays for the barotropic solution 
    5150#endif 
    5251 

  • trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r1241 r1438  
    11MODULE dynspg_ts 
    22   !!====================================================================== 
    3    !!                   ***  MODULE  dynspg_ts  *** 
    4    !! Ocean dynamics:  surface pressure gradient trend 
    5    !!====================================================================== 
    6    !! History :   9.0  !  04-12  (L. Bessieres, G. Madec)  Original code 
    7    !!             " "  !  05-11  (V. Garnier, G. Madec)  optimization 
    8    !!             9.0  !  06-08  (S. Masson)  distributed restart using iom 
    9    !!              "   !  08-01  (R. Benshila)  change averaging method 
    10    !!              "   !  07-07  (D. Storkey) calls to BDY routines 
    11    !!--------------------------------------------------------------------- 
     3   !! History :   1.0  !  04-12  (L. Bessieres, G. Madec)  Original code 
     4   !!              -   !  05-11  (V. Garnier, G. Madec)  optimization 
     5   !!              -   !  06-08  (S. Masson)  distributed restart using iom 
     6   !!             2.0  !  07-07  (D. Storkey) calls to BDY routines 
     7   !!              -   !  08-01  (R. Benshila)  change averaging method 
     8  !!--------------------------------------------------------------------- 
    129#if defined key_dynspg_ts   ||   defined key_esopa 
    1310   !!---------------------------------------------------------------------- 
    … …  
    1916   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file 
    2017   !!---------------------------------------------------------------------- 
    21    !! * Modules used 
    2218   USE oce             ! ocean dynamics and tracers 
    2319   USE dom_oce         ! ocean space and time domain 
    … …  
    4945   PUBLIC ts_rst      ! routine called by istate.F90 
    5046 
    51    REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne,   &  ! triad of coriolis parameter 
    52       &                             ftsw, ftse       ! (only used with een vorticity scheme) 
    53  
     47   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne   ! triad of coriolis parameter 
     48   REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme) 
    5449 
    5550   !! * Substitutions 
    5651#  include "domzgr_substitute.h90" 
    5752#  include "vectopt_loop_substitute.h90" 
    58    !!---------------------------------------------------------------------- 
    59    !!  OPA 9.0 , LOCEAN-IPSL (2005)  
     53   !!------------------------------------------------------------------------- 
     54   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    6055   !! $Id$ 
    61    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    62    !!---------------------------------------------------------------------- 
     56   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
     57   !!------------------------------------------------------------------------- 
    6358 
    6459CONTAINS 
    … …  
    8580      !!          (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b  
    8681      !!          and zvX_b (X specifying after, now or before). 
    87       !!      -3- Update of sea surface height from time averaged barotropic  
    88       !!          variables. 
    89       !!        - apply lateral boundary conditions on sshn. 
    90       !!      -4- The new general trend becomes : 
    91       !!          ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H 
     82      !!      -3- The new general trend becomes : 
     83      !!          ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H 
    9284      !! 
    9385      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 
    … …  
    9688      !!--------------------------------------------------------------------- 
    9789      INTEGER, INTENT( in )  ::   kt           ! ocean time-step index 
    98  
    99       !! * Local declarations 
     90      !! 
    10091      INTEGER  ::  ji, jj, jk, jit             ! dummy loop indices 
    10192      INTEGER  ::  icycle                      ! temporary scalar 
    … …  
    10798         zcu, zcv, zwx, zwy, zhdiv,         &  ! temporary arrays 
    10899         zua, zva, zub, zvb,                &  !     "        " 
    109          zssha_b, zua_b, zva_b,             &  !     "        " 
    110          zub_e, zvb_e,                      &  !     "        " 
    111          zun_e, zvn_e                          !     "        " 
     100         zua_e, zva_e, zssha_e,             &  !     "        " 
     101         zub_e, zvb_e, zsshb_e,             &  !     "        " 
     102         zu_sum, zv_sum, zssh_sum 
    112103      !! Variable volume 
    113104      REAL(wp), DIMENSION(jpi,jpj)     ::   & 
    114105         zspgu_1, zspgv_1, zsshun_e, zsshvn_e                     ! 2D workspace 
    115       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfse3un_e, zfse3vn_e  ! 3D workspace 
    116106      !!---------------------------------------------------------------------- 
    117107 
    118108      ! Arrays initialization 
    119109      ! --------------------- 
    120       zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0 
    121       zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0 
    122110      zhdiv(:,:) = 0.e0 
    123111 
    … …  
    133121         !                               ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 
    134122 
    135          ssha_e(:,:) = sshn(:,:) 
    136          ua_e(:,:)   = un_b(:,:) 
    137          va_e(:,:)   = vn_b(:,:) 
    138          hu_e(:,:)   = hu(:,:) 
    139          hv_e(:,:)   = hv(:,:) 
    140  
     123         zssha_e(:,:) = sshn(:,:) 
     124         zua_e  (:,:) = un_e(:,:) 
     125         zva_e  (:,:) = vn_e(:,:) 
     126         hu_e   (:,:) = hu  (:,:) 
     127         hv_e   (:,:) = hv  (:,:) 
    141128         IF( ln_dynvor_een ) THEN 
    142129            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0 
    … …  
    170157               zspgv_1(ji,jj) = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1) & 
    171158                  &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e2v(ji,jj) 
    172             END DO  
     159            END DO 
    173160         END DO 
    174161 
    … …  
    193180      zfact2 = 0.5 * 0.5 
    194181      zraur  = 1. / rauw                                    ! 1 / volumic mass of pure water 
    195        
     182 
    196183      ! ----------------------------------------------------------------------------- 
    197184      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step) 
    … …  
    215202               zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk) 
    216203               !                                                           ! Vertically integrated transports (before) 
    217                zub(ji,1) = zub(ji,1) + fse3u(ji,1,jk) * ub(ji,1,jk) 
    218                zvb(ji,1) = zvb(ji,1) + fse3v(ji,1,jk) * vb(ji,1,jk) 
     204               zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 
     205               zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 
    219206               !                                                           ! Planetary vorticity transport fluxes (now) 
    220207               zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk) 
    … …  
    228215            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    229216            !                                                           ! Vertically integrated transports (before) 
    230             zub(:,:) = zub(:,:) + fse3u(:,:,jk) * ub(:,:,jk) 
    231             zvb(:,:) = zvb(:,:) + fse3v(:,:,jk) * vb(:,:,jk) 
     217            zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 
     218            zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 
    232219            !                                                           ! Planetary vorticity (now) 
    233220            zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    … …  
    304291      !---------------- 
    305292      ! Number of iteration of the barotropic loop 
    306       icycle = 3  * nn_baro / 2 
     293      icycle = 2  * nn_baro + 1 
    307294 
    308295      ! variables for the barotropic equations 
    309       sshb_e (:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now) 
    310       sshn_e (:,:) = sshn_b(:,:) 
    311       zub_e  (:,:) = un_b  (:,:)       ! barotropic transports issued from the barotropic equations (before and now) 
    312       zvb_e  (:,:) = vn_b  (:,:) 
    313       zun_e  (:,:) = un_b  (:,:) 
    314       zvn_e  (:,:) = vn_b  (:,:) 
    315       zssha_b(:,:) = 0.e0 
    316       zua_b  (:,:) = 0.e0 
    317       zva_b  (:,:) = 0.e0 
    318       hu_e   (:,:) = hu   (:,:)     ! (barotropic) ocean depth at u-point 
    319       hv_e   (:,:) = hv   (:,:)     ! (barotropic) ocean depth at v-point 
    320       IF( lk_vvl ) THEN 
    321          zsshun_e(:,:)    = sshu (:,:)     ! (barotropic) sea surface height (now) at u-point 
    322          zsshvn_e(:,:)    = sshv (:,:)     ! (barotropic) sea surface height (now) at v-point 
    323          zfse3un_e(:,:,:) = fse3u(:,:,:)   ! (barotropic) scale factors  at u-point 
    324          zfse3un_e(:,:,:) = fse3v(:,:,:)   ! (barotropic) scale factors  at v-point 
     296      zu_sum  (:,:) = 0.e0 
     297      zv_sum  (:,:) = 0.e0 
     298      zssh_sum(:,:) = 0.e0 
     299      hu_e    (:,:) = hu    (:,:)      ! (barotropic) ocean depth at u-point 
     300      hv_e    (:,:) = hv    (:,:)      ! (barotropic) ocean depth at v-point 
     301      zsshb_e (:,:) = sshn_e(:,:)      ! (barotropic) sea surface height (before and now) 
     302      ! vertical sum 
     303      un_e  (:,:) = 0.e0 
     304      vn_e  (:,:) = 0.e0 
     305      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
     306         DO jk = 1, jpkm1 
     307            DO ji = 1, jpij 
     308               un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 
     309               vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 
     310            END DO 
     311         END DO 
     312      ELSE                             ! No  vector opt. 
     313         DO jk = 1, jpkm1 
     314            un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
     315            vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
     316         END DO 
    325317      ENDIF 
     318      zub_e (:,:) = un_e(:,:) 
     319      zvb_e (:,:) = vn_e(:,:) 
     320 
    326321 
    327322      ! set ssh corrections to 0 
    … …  
    352347         DO jj = 2, jpjm1 
    353348            DO ji = fs_2, fs_jpim1   ! vector opt. 
    354                zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * zun_e(ji  ,jj)              & 
    355                   &            -e2u(ji-1,jj  ) * zun_e(ji-1,jj)              & 
    356                   &            +e1v(ji  ,jj  ) * zvn_e(ji  ,jj)              & 
    357                   &            -e1v(ji  ,jj-1) * zvn_e(ji  ,jj-1) )          & 
     349               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * un_e(ji  ,jj)              & 
     350                  &            -e2u(ji-1,jj  ) * un_e(ji-1,jj)              & 
     351                  &            +e1v(ji  ,jj  ) * vn_e(ji  ,jj)              & 
     352                  &            -e1v(ji  ,jj-1) * vn_e(ji  ,jj-1) )          & 
    358353                  &           / (e1t(ji,jj)*e2t(ji,jj)) 
    359354            END DO 
    … …  
    370365 
    371366#if defined key_bdy 
    372             DO jj = 1, jpj 
    373                DO ji = 1, jpi 
    374                   zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 
    375                END DO 
    376             END DO 
     367         DO jj = 1, jpj 
     368            DO ji = 1, jpi 
     369               zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 
     370            END DO 
     371         END DO 
    377372#endif 
    378373 
    … …  
    381376         DO jj = 2, jpjm1 
    382377            DO ji = fs_2, fs_jpim1   ! vector opt. 
    383                ssha_e(ji,jj) = ( sshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  & 
    384             &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
     378               zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  & 
     379                  &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
    385380            END DO 
    386381         END DO 
    … …  
    388383         ! evolution of the barotropic transport ( following the vorticity scheme used) 
    389384         ! ---------------------------------------------------------------------------- 
    390          zwx(:,:) = e2u(:,:) * zun_e(:,:) 
    391          zwy(:,:) = e1v(:,:) * zvn_e(:,:) 
     385         zwx(:,:) = e2u(:,:) * un_e(:,:) 
     386         zwy(:,:) = e1v(:,:) * vn_e(:,:) 
    392387 
    393388         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
    … …  
    396391                  ! surface pressure gradient 
    397392                  IF( lk_vvl) THEN 
    398                      zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 
    399                         &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 
    400                      zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 
    401                         &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 
     393                     zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) & 
     394                        &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj) 
     395                     zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) & 
     396                        &            - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj) 
    402397                  ELSE 
    403398                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 
    … …  
    412407                  zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) 
    413408                  ! after transports 
    414                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    415                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     409                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
     410                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
    416411               END DO 
    417412            END DO 
    … …  
    422417                  ! surface pressure gradient 
    423418                  IF( lk_vvl) THEN 
    424                      zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 
    425                         &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 
    426                      zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 
    427                         &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 
     419                     zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) & 
     420                        &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj) 
     421                     zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) & 
     422                        &            - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj) 
    428423                  ELSE 
    429424                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 
    … …  
    432427                  ! enstrophy conserving formulation for planetary vorticity term 
    433428                  zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   & 
    434                                  + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     429                     + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    435430                  zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   & 
    436                                  + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     431                     + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    437432                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) 
    438433                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) 
    439434                  ! after transports 
    440                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    441                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     435                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
     436                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
    442437               END DO 
    443438            END DO 
    … …  
    449444                  ! surface pressure gradient 
    450445                  IF( lk_vvl) THEN 
    451                      zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 
    452                         &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 
    453                      zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 
    454                         &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 
     446                     zspgu = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  ) & 
     447                        &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj) 
     448                     zspgv = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) & 
     449                        &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj) 
    455450                  ELSE 
    456451                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 
    … …  
    463458                     &                             + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    464459                  ! after transports 
    465                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    466                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     460                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
     461                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
    467462               END DO 
    468463            END DO 
    … …  
    477472 
    478473         ! ... Boundary conditions on ua_e, va_e, ssha_e 
    479          CALL lbc_lnk( ua_e  , 'U', -1. ) 
    480          CALL lbc_lnk( va_e  , 'V', -1. ) 
    481          CALL lbc_lnk( ssha_e, 'T',  1. ) 
     474         CALL lbc_lnk( zua_e  , 'U', -1. ) 
     475         CALL lbc_lnk( zva_e  , 'V', -1. ) 
     476         CALL lbc_lnk( zssha_e, 'T',  1. ) 
    482477 
    483478         ! temporal sum 
    484479         !------------- 
    485          IF( jit >= nn_baro / 2 ) THEN 
    486             zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:) 
    487             zua_b  (:,:) = zua_b  (:,:) + ua_e  (:,:) 
    488             zva_b  (:,:) = zva_b  (:,:) + va_e  (:,:)  
    489          ENDIF 
     480         zu_sum  (:,:) = zu_sum  (:,:) + zua_e  (:,:) 
     481         zv_sum  (:,:) = zv_sum  (:,:) + zva_e  (:,:)  
     482         zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:)  
    490483 
    491484         ! Time filter and swap of dynamics arrays 
    492485         ! --------------------------------------- 
    493486         IF( jit == 1 ) THEN   ! Euler (forward) time stepping 
    494             sshb_e (:,:) = sshn_e(:,:) 
    495             zub_e  (:,:) = zun_e (:,:) 
    496             zvb_e  (:,:) = zvn_e (:,:) 
    497             sshn_e (:,:) = ssha_e(:,:) 
    498             zun_e  (:,:) = ua_e  (:,:) 
    499             zvn_e  (:,:) = va_e  (:,:) 
     487            zsshb_e(:,:) = sshn_e (:,:) 
     488            zub_e  (:,:) = un_e  (:,:) 
     489            zvb_e  (:,:) = vn_e  (:,:) 
     490            sshn_e (:,:) = zssha_e(:,:) 
     491            un_e   (:,:) = zua_e  (:,:) 
     492            vn_e   (:,:) = zva_e  (:,:) 
    500493         ELSE                                        ! Asselin filtering 
    501             sshb_e (:,:) = atfp * ( sshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 
    502             zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e  (:,:) 
    503             zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e  (:,:) 
    504             sshn_e (:,:) = ssha_e(:,:) 
    505             zun_e  (:,:) = ua_e  (:,:) 
    506             zvn_e  (:,:) = va_e  (:,:) 
     494            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:) 
     495            zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * un_e  (:,:) 
     496            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * vn_e  (:,:) 
     497            sshn_e (:,:) = zssha_e(:,:) 
     498            un_e   (:,:) = zua_e  (:,:) 
     499            vn_e   (:,:) = zva_e  (:,:) 
    507500         ENDIF 
    508501 
    509          IF( lk_vvl ) THEN ! Variable volume 
     502         IF( lk_vvl ) THEN ! Variable volume   ! needed for BDY ??? 
    510503 
    511504            ! Sea surface elevation 
    … …  
    528521 
    529522            ! Boundaries conditions 
    530             CALL lbc_lnk( zsshun_e, 'U', 1. )    ;    CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
    531  
    532             ! Scale factors at before and after time step 
     523            CALL lbc_lnk( zsshun_e, 'U', 1. )    ;  CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
     524 
     525            ! Ocean depth at U- and V-points 
    533526            ! ------------------------------------------- 
    534             CALL dom_vvl_sf( zsshun_e, 'U', zfse3un_e ) ;   CALL dom_vvl_sf( zsshvn_e, 'V', zfse3vn_e ) 
    535  
    536             ! Ocean depth at U- and V-points 
    537             hu_e(:,:) = 0.e0 
    538             hv_e(:,:) = 0.e0 
    539  
    540             DO jk = 1, jpk 
    541                hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk) 
    542                hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk) 
    543             END DO 
    544  
    545          ENDIF ! End variable volume case 
    546  
     527            hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:) 
     528            hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:) 
     529 
     530            ! 
     531         ENDIF 
    547532         !                                                 ! ==================== ! 
    548533      END DO                                               !        end loop      ! 
    549534      !                                                    ! ==================== ! 
    550535 
    551  
    552536      ! Time average of after barotropic variables 
    553       zcoef =  1.e0 / ( nn_baro + 1 )  
    554       zssha_b(:,:) = zcoef * zssha_b(:,:)  
    555       zua_b  (:,:) = zcoef *  zua_b (:,:)  
    556       zva_b  (:,:) = zcoef *  zva_b (:,:)  
     537      zcoef =  1.e0 / ( 2 * nn_baro  + 1 )  
     538      un_e  (:,:) = zcoef *  zu_sum  (:,:)  
     539      vn_e  (:,:) = zcoef *  zv_sum  (:,:)  
     540      sshn_e(:,:) = zcoef *  zssh_sum(:,:)  
     541 
    557542#if defined key_obc 
    558543      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) 
    … …  
    561546      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:) 
    562547#endif 
    563       
    564  
    565       ! --------------------------------------------------------------------------- 
    566       ! Phase 3 : Update sea surface height from time averaged barotropic variables 
    567       ! --------------------------------------------------------------------------- 
    568  
    569   
    570       ! Horizontal divergence of time averaged barotropic transports 
    571       !------------------------------------------------------------- 
    572       IF( .NOT. lk_vvl ) THEN 
    573          DO jj = 2, jpjm1 
    574             DO ji = fs_2, fs_jpim1   ! vector opt. 
    575                zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     & 
    576                   &            +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   & 
    577                   &           / ( e1t(ji,jj) * e2t(ji,jj) ) 
    578             END DO 
    579          END DO 
    580       ENDIF 
    581  
    582 #if defined key_obc && ! defined key_vvl 
    583       ! open boundaries (div must be zero behind the open boundary) 
    584       !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    585       IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east 
    586       IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west 
    587       IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north 
    588       IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south 
    589 #endif 
    590  
    591 #if defined key_bdy 
    592          DO jj = 1, jpj 
    593            DO ji = 1, jpi 
    594              zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 
    595            END DO 
    596          END DO 
    597 #endif 
    598  
    599       ! sea surface height 
    600       !------------------- 
    601       IF( lk_vvl ) THEN 
    602          sshbb(:,:) = sshb(:,:) 
    603          sshb (:,:) = sshn(:,:) 
    604          sshn (:,:) = ssha(:,:) 
    605       ELSE 
    606          sshb (:,:) = sshn(:,:) 
    607          sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    608       ENDIF 
    609  
    610       ! ... Boundary conditions on sshn 
    611       IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. ) 
    612  
    613548 
    614549      ! ----------------------------------------------------------------------------- 
    615       ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step) 
     550      ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step) 
    616551      ! ----------------------------------------------------------------------------- 
    617552 
    618       ! Swap on time averaged barotropic variables 
    619       ! ------------------------------------------ 
    620       sshb_b(:,:) = sshn_b (:,:) 
    621       IF ( neuler == 0 .AND. kt == nit000 ) zssha_b(:,:) = sshn(:,:) 
    622       sshn_b(:,:) = zssha_b(:,:) 
    623       un_b  (:,:) = zua_b  (:,:)  
    624       vn_b  (:,:) = zva_b  (:,:)  
    625     
     553 
     554 
    626555      ! add time averaged barotropic coriolis and surface pressure gradient 
    627556      ! terms to the general momentum trend 
    628557      ! -------------------------------------------------------------------- 
    629558      DO jk=1,jpkm1 
    630          ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b 
    631          va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b 
     559         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b 
     560         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b 
    632561      END DO 
    633562 
    … …  
    636565      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' ) 
    637566 
    638       ! print sum trends (used for debugging) 
    639       IF( ln_ctl )     CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask ) 
    640567      ! 
    641568   END SUBROUTINE dyn_spg_ts 
    … …  
    655582      ! 
    656583      IF( TRIM(cdrw) == 'READ' ) THEN 
    657          IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
    658             CALL iom_get( numror, jpdom_autoglo, 'sshb'  , sshb(:,:)   ) 
    659             CALL iom_get( numror, jpdom_autoglo, 'sshn'  , sshn(:,:)   ) 
    660             IF( neuler == 0 ) sshb(:,:) = sshn(:,:) 
     584         IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN 
     585            CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) )   ! free surface and 
     586            CALL iom_get( numror, jpdom_autoglo, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued 
     587            CALL iom_get( numror, jpdom_autoglo, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop 
    661588         ELSE 
    662             IF( nn_rstssh == 1 ) THEN   
    663                sshb(:,:) = 0.e0 
    664                sshn(:,:) = 0.e0 
    665             ENDIF 
    666          ENDIF 
    667          IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 
    668             CALL iom_get( numror, jpdom_autoglo, 'sshb_b', sshb_b(:,:) )   ! free surface issued 
    669             CALL iom_get( numror, jpdom_autoglo, 'sshn_b', sshn_b(:,:) )   ! from time-splitting loop 
    670             CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued 
    671             CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
    672             IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:) 
    673          ELSE 
    674             sshb_b(:,:) = sshb(:,:) 
    675             sshn_b(:,:) = sshn(:,:) 
    676             un_b  (:,:) = 0.e0 
    677             vn_b  (:,:) = 0.e0 
     589            sshn_e(:,:) = sshn(:,:) 
     590            un_e  (:,:) = 0.e0 
     591            vn_e  (:,:) = 0.e0 
    678592            ! vertical sum 
    679593            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    680594               DO jk = 1, jpkm1 
    681595                  DO ji = 1, jpij 
    682                      un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 
    683                      vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 
     596                     un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 
     597                     vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 
    684598                  END DO 
    685599               END DO 
    686600            ELSE                             ! No  vector opt. 
    687601               DO jk = 1, jpkm1 
    688                   un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    689                   vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
     602                  un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
     603                  vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    690604               END DO 
    691605            ENDIF 
    692606         ENDIF 
    693607      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    694          CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) ) 
    695          CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) ) 
    696          CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) )   ! free surface issued 
    697          CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) )   ! from barotropic loop 
    698          CALL iom_rstput( kt, nitrst, numrow, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued 
    699          CALL iom_rstput( kt, nitrst, numrow, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
     608         CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) )   ! free surface and 
     609         CALL iom_rstput( kt, nitrst, numrow, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued 
     610         CALL iom_rstput( kt, nitrst, numrow, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop 
    700611      ENDIF 
    701612      ! 

  • trunk/NEMO/OPA_SRC/DYN/dynvor.F90

    r1152 r1438  
    55   !!                 planetary vorticity trends 
    66   !!====================================================================== 
    7    !! History :  1.0  !  89-12  (P. Andrich)  vor_ens: Original code 
    8    !!            5.0  !  91-11  (G. Madec) vor_ene, vor_mix: Original code 
    9    !!            6.0  !  96-01  (G. Madec)  s-coord, suppress work arrays 
    10    !!            8.5  !  02-08  (G. Madec)  F90: Free form and module 
    11    !!            8.5  !  04-02  (G. Madec)  vor_een: Original code 
    12    !!            9.0  !  03-08  (G. Madec)  vor_ctl: Original code 
    13    !!            9.0  !  05-11  (G. Madec)  dyn_vor: Original code (new step architecture) 
    14    !!            9.0  !  06-11  (G. Madec)  flux form advection: add metric term 
     7   !! History :  OPA  !  1989-12  (P. Andrich)  vor_ens: Original code 
     8   !!            5.0  !  1991-11  (G. Madec) vor_ene, vor_mix: Original code 
     9   !!            6.0  !  1996-01  (G. Madec)  s-coord, suppress work arrays 
     10   !!            8.5  !  2002-08  (G. Madec)  F90: Free form and module 
     11   !!   NEMO     1.0  !  2004-02  (G. Madec)  vor_een: Original code 
     12   !!             -   !  2003-08  (G. Madec)  add vor_ctl 
     13   !!             -   !  2005-11  (G. Madec)  add dyn_vor (new step architecture) 
     14   !!            2.0  !  2006-11  (G. Madec)  flux form advection: add metric term 
     15   !!            3.2  !  2009-04  (R. Benshila)  vvl: correction of een scheme 
    1516   !!---------------------------------------------------------------------- 
    1617 
    … …  
    3738   PUBLIC   dyn_vor   ! routine called by step.F90 
    3839 
    39    !!* Namelist nam_dynvor: vorticity term 
     40   !                                             !!* Namelist nam_dynvor: vorticity term 
    4041   LOGICAL, PUBLIC ::   ln_dynvor_ene = .FALSE.   !: energy conserving scheme 
    4142   LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme 
    … …  
    5253#  include "vectopt_loop_substitute.h90" 
    5354   !!---------------------------------------------------------------------- 
    54    !!   OPA 9.0 , LOCEAN-IPSL (2006)  
     55   !! NEMO/OPA 3,2 , LOCEAN-IPSL (2009)  
    5556   !! $Id$ 
    5657   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    174175      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, & 
    175176         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    176  
     177      ! 
    177178   END SUBROUTINE dyn_vor 
    178179 
    … …  
    206207      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    207208      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    208          !                                                        ! =nrvm (relative vorticity or metric) 
     209      !                                                           ! =nrvm (relative vorticity or metric) 
    209210      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    210211      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    … …  
    222223      ENDIF 
    223224 
    224       ! Local constant initialization 
    225       zfact2 = 0.5 * 0.5 
     225      zfact2 = 0.5 * 0.5      ! Local constant initialization 
    226226 
    227227!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
    … …  
    229229      DO jk = 1, jpkm1                                 ! Horizontal slab 
    230230         !                                             ! =============== 
     231         ! 
    231232         ! Potential vorticity and horizontal fluxes 
    232233         ! ----------------------------------------- 
    … …  
    315316      INTEGER, INTENT(in) ::   kt   ! ocean timestep index 
    316317      !! 
    317       INTEGER ::   ji, jj, jk   ! dummy loop indices 
     318      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    318319      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! temporary scalars 
    319320      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !    "         " 
    … …  
    327328      ENDIF 
    328329 
    329       ! Local constant initialization 
    330       zfact1 = 0.5 * 0.25 
     330      zfact1 = 0.5 * 0.25      ! Local constant initialization 
    331331      zfact2 = 0.5 * 0.5 
    332332 
    … …  
    335335      DO jk = 1, jpkm1                                 ! Horizontal slab 
    336336         !                                             ! =============== 
    337  
     337         ! 
    338338         ! Relative and planetary potential vorticity and horizontal fluxes 
    339339         ! ---------------------------------------------------------------- 
    … …  
    438438      ENDIF 
    439439 
    440       ! Local constant initialization 
    441       zfact1 = 0.5 * 0.25 
     440      zfact1 = 0.5 * 0.25      ! Local constant initialization 
    442441 
    443442!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
    … …  
    445444      DO jk = 1, jpkm1                                 ! Horizontal slab 
    446445         !                                             ! =============== 
     446         ! 
    447447         ! Potential vorticity and horizontal fluxes 
    448448         ! ----------------------------------------- 
    … …  
    465465                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    466466                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    467                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               & 
     467                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                & 
    468468                       &       ) 
    469469               END DO 
    470470            END DO 
    471471         END SELECT 
    472  
     472         ! 
    473473         IF( ln_sco ) THEN 
    474474            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop  
    … …  
    487487            END DO 
    488488         ENDIF 
    489  
     489         ! 
    490490         ! Compute and add the vorticity term trend 
    491491         ! ---------------------------------------- 
    … …  
    514514      !! 
    515515      !! ** Method  :   Trend evaluated using now fields (centered in time)  
    516       !!      and the Arakawa and Lamb (19XX) flux form formulation : conserves  
     516      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves  
    517517      !!      both the horizontal kinetic energy and the potential enstrophy 
    518       !!      when horizontal divergence is zero. 
    519       !!      The trend of the vorticity term is given by: 
    520       !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 
    521       !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 
    522       !!      Add this trend to the general momentum trend (ua,va): 
    523       !!          (ua,va) = (ua,va) + ( voru , vorv ) 
     518      !!      when horizontal divergence is zero (see the NEMO documentation) 
     519      !!      Add this trend to the general momentum trend (ua,va). 
    524520      !! 
    525521      !! ** Action : - Update (ua,va) with the now vorticity term trend 
    … …  
    531527      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    532528      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    533          !                                                        ! =nrvm (relative vorticity or metric) 
     529      !                                                           ! =nrvm (relative vorticity or metric) 
    534530      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    535531      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    536532      !! 
    537       INTEGER ::   ji, jj, jk          ! dummy loop indices 
     533      INTEGER  ::   ji, jj, jk         ! dummy loop indices 
    538534      REAL(wp) ::   zfac12, zua, zva   ! temporary scalars 
    539535      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz            ! temporary 2D workspace 
    540536      REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse   ! temporary 3D workspace 
     537#if defined key_vvl 
     538      REAL(wp), DIMENSION(jpi,jpj,jpk)       ::   ze3f 
     539#else 
    541540      REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE ::   ze3f 
     541#endif 
    542542      !!---------------------------------------------------------------------- 
    543543 
    … …  
    546546         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 
    547547         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    548  
     548      ENDIF 
     549 
     550      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t) 
    549551         DO jk = 1, jpk 
    550552            DO jj = 1, jpjm1 
    … …  
    559561      ENDIF 
    560562 
    561       ! Local constant initialization 
    562       zfac12 = 1.e0 / 12.e0 
     563      zfac12 = 1.e0 / 12.e0      ! Local constant initialization 
    563564 
    564565       
    … …  
    571572         ! ----------------------------------------- 
    572573         SELECT CASE( kvor )      ! vorticity considered 
    573          CASE ( 1 )   ;   zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)   ! planetary vorticity (Coriolis) 
    574          CASE ( 2 )   ;   zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)   ! relative  vorticity 
     574         CASE ( 1 )                                                ! planetary vorticity (Coriolis) 
     575            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk) 
     576         CASE ( 2 )                                                ! relative  vorticity 
     577            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 
    575578         CASE ( 3 )                                                ! metric term 
    576579            DO jj = 1, jpjm1 
    … …  
    581584               END DO 
    582585            END DO 
    583          CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) ! total (relative + planetary vorticity) 
     586         CASE ( 4 )                                                ! total (relative + planetary vorticity) 
     587            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 
    584588         CASE ( 5 )                                                ! total (coriolis + metric) 
    585589            DO jj = 1, jpjm1 
    … …  
    588592                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    589593                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    590                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               & 
     594                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                & 
    591595                       &       ) * ze3f(ji,jj,jk) 
    592596               END DO 
    … …  
    599603         ! Compute and add the vorticity term trend 
    600604         ! ---------------------------------------- 
    601          jj=2 
    602          ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 
     605         jj = 2 
     606         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;  ztsw(1,:) = 0 
    603607         DO ji = 2, jpi    
    604608               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
    … …  
    636640      !! 
    637641      !! ** Purpose :   Control the consistency between cpp options for 
    638       !!      tracer advection schemes 
     642      !!              tracer advection schemes 
    639643      !!---------------------------------------------------------------------- 
    640644      INTEGER ::   ioptio          ! temporary integer 

  • trunk/NEMO/OPA_SRC/DYN/wzvmod.F90

    r1241 r1438  
    11MODULE wzvmod 
     2!! MODULE sshwzv    
    23   !!============================================================================== 
    3    !!                       ***  MODULE  wzvmod  *** 
    4    !! Ocean diagnostic variable : vertical velocity 
     4   !!                       ***  MODULE  sshwzv  *** 
     5   !! Ocean dynamics : sea surface height and vertical velocity 
    56   !!============================================================================== 
    6    !! History :   5.0  !  90-10  (C. Levy, G. Madec)  Original code 
    7    !!             7.0  !  96-01  (G. Madec)  Statement function for e3 
    8    !!             8.5  !  02-07  (G. Madec)  Free form, F90 
    9    !!              "   !  07-07  (D. Storkey) Zero zhdiv at open boundary (BDY)  
    10    !!---------------------------------------------------------------------- 
    11    !!   wzv        : Compute the vertical velocity 
    12    !!---------------------------------------------------------------------- 
    13    !! * Modules used 
     7   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code 
     8   !!---------------------------------------------------------------------- 
     9 
     10   !!---------------------------------------------------------------------- 
     11   !!   ssh_wzv        : after ssh & now vertical velocity 
     12   !!   ssh_nxt        : filter ans swap the ssh arrays 
     13   !!   ssh_rst        : read/write ssh restart fields in the ocean restart file 
     14   !!---------------------------------------------------------------------- 
    1415   USE oce             ! ocean dynamics and tracers variables 
    1516   USE dom_oce         ! ocean space and time domain variables  
    1617   USE sbc_oce         ! surface boundary condition: ocean 
    1718   USE domvvl          ! Variable volume 
     19   USE iom             ! I/O library 
     20   USE restart         ! only for lrst_oce 
    1821   USE in_out_manager  ! I/O manager 
    1922   USE prtctl          ! Print control 
    2023   USE phycst 
    21    USE bdy_oce         ! unstructured open boundaries 
    2224   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    2325   USE obc_par         ! open boundary cond. parameter 
    … …  
    2729   PRIVATE 
    2830 
    29    !! * Routine accessibility 
    30    PUBLIC wzv          ! routine called by step.F90 and inidtr.F90 
     31   PUBLIC   ssh_wzv    ! called by step.F90 
     32   PUBLIC   ssh_nxt    ! called by step.F90 
    3133 
    3234   !! * Substitutions 
    3335#  include "domzgr_substitute.h90" 
    34    !!---------------------------------------------------------------------- 
    35    !!  OPA 9.0 , LOCEAN-IPSL (2005)  
     36#  include "vectopt_loop_substitute.h90" 
     37 
     38   !!---------------------------------------------------------------------- 
     39   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    3640   !! $Id$ 
    3741   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    4044CONTAINS 
    4145 
    42    SUBROUTINE wzv( kt ) 
    43       !!---------------------------------------------------------------------- 
    44       !!                    ***  ROUTINE wzv  *** 
    45       !! 
    46       !! ** Purpose :   Compute the now vertical velocity after the array swap 
    47       !! 
    48       !! ** Method  :   Using the incompressibility hypothesis, the vertical 
    49       !!      velocity is computed by integrating the horizontal divergence  
    50       !!      from the bottom to the surface. 
    51       !!        The boundary conditions are w=0 at the bottom (no flux) and, 
    52       !!      in regid-lid case, w=0 at the sea surface. 
    53       !! 
    54       !! ** action  :   wn array : the now vertical velocity 
    55       !!---------------------------------------------------------------------- 
    56       !! * Arguments 
    57       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    58  
    59       !! * Local declarations 
    60       INTEGER  ::           jk           ! dummy loop indices 
    61       !! Variable volume 
    62       INTEGER  ::   ji, jj               ! dummy loop indices 
    63       REAL(wp) ::   z2dt, zraur          ! temporary scalar 
    64       REAL(wp), DIMENSION (jpi,jpj) ::   zssha, zun, zvn, zhdiv 
    65 #if defined key_bdy 
    66       INTEGER  ::     jgrd, jb           ! temporary scalars 
    67 #endif 
     46   SUBROUTINE ssh_wzv( kt )  
     47      !!---------------------------------------------------------------------- 
     48      !!                ***  ROUTINE ssh_wzv  *** 
     49      !!                    
     50      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity 
     51      !!              and update the now vertical coordinate (lk_vvl=T). 
     52      !! 
     53      !! ** Method  : - 
     54      !! 
     55      !!              - Using the incompressibility hypothesis, the vertical  
     56      !!      velocity is computed by integrating the horizontal divergence   
     57      !!      from the bottom to the surface minus the scale factor evolution. 
     58      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     59      !! 
     60      !! ** action  :   ssha    : after sea surface height 
     61      !!                wn      : now vertical velocity 
     62      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height 
     63      !!                          at u-, v-, f-point s 
     64      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 
     65      !!---------------------------------------------------------------------- 
     66      INTEGER, INTENT(in) ::   kt   ! time step 
     67      !! 
     68      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     69      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars 
     70      REAL(wp) ::   z2dt, zraur     ! temporary scalars 
     71      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace 
    6872      !!---------------------------------------------------------------------- 
    6973 
    7074      IF( kt == nit000 ) THEN 
    7175         IF(lwp) WRITE(numout,*) 
    72          IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.' 
    73          IF(lwp) WRITE(numout,*) '~~~~~~~ '  
    74  
    75          ! bottom boundary condition: w=0 (set once for all) 
    76          wn(:,:,jpk) = 0.e0 
    77       ENDIF 
    78  
    79       IF( lk_vvl ) THEN                ! Variable volume 
    80          ! 
    81          z2dt = 2. * rdt                                       ! time step: leap-frog 
    82          IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest 
    83          zraur  = 1. / rauw 
    84  
    85          ! Vertically integrated quantities 
    86          ! -------------------------------- 
    87          zun(:,:) = 0.e0 
    88          zvn(:,:) = 0.e0 
    89          ! 
    90          DO jk = 1, jpkm1             ! Vertically integrated transports (now) 
    91             zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    92             zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    93          END DO 
    94  
    95          ! Horizontal divergence of barotropic transports 
    96          !-------------------------------------------------- 
    97          zhdiv(:,:) = 0.e0 
    98          DO jj = 2, jpjm1 
    99             DO ji = 2, jpim1   ! vector opt. 
    100                zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )     & 
    101                   &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )     & 
    102                   &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )     & 
    103                   &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )   & 
    104                   &           / ( e1t(ji,jj) * e2t(ji,jj) ) 
     76         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' 
     77         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     78         ! 
     79         CALL ssh_rst( nit000, 'READ' )       ! read or initialize sshb and sshn 
     80         ! 
     81         wn(:,:,jpk) = 0.e0                   ! bottom boundary condition: w=0 (set once for all) 
     82         ! 
     83         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only) 
     84            DO jj = 1, jpjm1 
     85               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible 
     86                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     87                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     88                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
     89                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
     90                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
     91                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
     92                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
     93                  sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 & 
     94                     &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) 
     95                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
     96                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
     97                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
     98                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
     99                  sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 & 
     100                     &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 
     101               END DO 
     102            END DO 
     103            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
     104            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
     105            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. ) 
     106         ENDIF 
     107         ! 
     108      ENDIF 
     109 
     110      ! set time step size (Euler/Leapfrog) 
     111      z2dt = 2. * rdt  
     112      IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 
     113 
     114      zraur = 1. / rauw 
     115 
     116      !                                           !------------------------------! 
     117      !                                           !   After Sea Surface Height   ! 
     118      !                                           !------------------------------! 
     119      zhdiv(:,:) = 0.e0 
     120      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
     121        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 
     122      END DO 
     123 
     124      !                                                ! Sea surface elevation time stepping 
     125      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     126 
     127#if defined key_obc 
     128# if defined key_agrif 
     129      IF ( Agrif_Root() ) THEN  
     130# endif 
     131         ssha(:,:) = ssha(:,:) * obctmsk(:,:) 
     132         CALL lbc_lnk( ssha, 'T', 1. )  ! absolutly compulsory !! (jmm) 
     133# if defined key_agrif 
     134      ENDIF 
     135# endif 
     136#endif 
     137 
     138      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only) 
     139      IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
     140         DO jj = 1, jpjm1 
     141            DO ji = 1, fs_jpim1   ! Vector Opt. 
     142               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
     143                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     & 
     144                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
     145               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
     146                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
     147                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     148               sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)   &   ! Caution : fmask not used 
     149                  &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 & 
     150                  &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) ) 
    105151            END DO 
    106152         END DO 
    107  
    108 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 
    109          ! open boundaries (div must be zero behind the open boundary) 
    110          !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    111          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east 
    112          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west 
    113          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north 
    114          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south 
    115 #endif 
    116  
    117 #if defined key_bdy 
    118          jgrd=1 !: tracer grid. 
    119          DO jb = 1, nblenrim(jgrd) 
    120            ji = nbi(jb,jgrd) 
    121            jj = nbj(jb,jgrd) 
    122            zhdiv(ji,jj) = 0.e0 
     153         CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions 
     154         CALL lbc_lnk( sshv_a, 'V', 1. ) 
     155         CALL lbc_lnk( sshf_a, 'F', 1. ) 
     156      ENDIF 
     157 
     158      !                                           !------------------------------! 
     159      !                                           !     Now Vertical Velocity    ! 
     160      !                                           !------------------------------! 
     161      !                                                ! integrate from the bottom the hor. divergence 
     162      DO jk = jpkm1, 1, -1 
     163         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   & 
     164              &                    - (  fse3t_a(:,:,jk)                   & 
     165              &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 
     166      END DO 
     167 
     168      !                                           !------------------------------! 
     169      !                                           !  Update Now Vertical coord.  ! 
     170      !                                           !------------------------------! 
     171      IF( lk_vvl ) THEN                           ! only in vvl case) 
     172         !                                             ! now local depth and scale factors (stored in fse3. arrays) 
     173         DO jk = 1, jpkm1 
     174            fsdept(:,:,jk) = fsdept_n(:,:,jk)               ! depths 
     175            fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 
     176            fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 
     177            ! 
     178            fse3t (:,:,jk) = fse3t_n (:,:,jk)               ! vertical scale factors 
     179            fse3u (:,:,jk) = fse3u_n (:,:,jk) 
     180            fse3v (:,:,jk) = fse3v_n (:,:,jk) 
     181            fse3f (:,:,jk) = fse3f_n (:,:,jk) 
     182            fse3w (:,:,jk) = fse3w_n (:,:,jk) 
     183            fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 
     184            fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 
    123185         END DO 
    124 #endif 
    125  
    126          CALL lbc_lnk( zhdiv, 'T', 1. ) 
    127  
    128          ! Sea surface elevation time stepping 
    129          ! ----------------------------------- 
    130          zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 
    131  
    132          ! Vertical velocity computed from bottom 
    133          ! -------------------------------------- 
    134          DO jk = jpkm1, 1, -1 
    135             wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) & 
    136               &        - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt 
    137          END DO 
    138  
    139       ELSE                             ! Fixed volume  
    140  
    141          ! Vertical velocity computed from bottom 
    142          ! -------------------------------------- 
    143          DO jk = jpkm1, 1, -1 
    144             wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) 
    145          END DO 
    146        
    147       ENDIF  
    148  
    149       IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn) 
    150  
    151    END SUBROUTINE wzv 
     186         !                                             ! ocean depth (at u- and v-points) 
     187         hu(:,:) = hu_0(:,:) + sshu_n(:,:) 
     188         hv(:,:) = hv_0(:,:) + sshv_n(:,:) 
     189         !                                             ! masked inverse of the ocean depth (at u- and v-points) 
     190         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 
     191         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 
     192 
     193      ENDIF 
     194      ! 
     195   END SUBROUTINE ssh_wzv 
     196 
     197 
     198   SUBROUTINE ssh_nxt( kt ) 
     199      !!---------------------------------------------------------------------- 
     200      !!                    ***  ROUTINE ssh_nxt  *** 
     201      !! 
     202      !! ** Purpose :   achieve the sea surface  height time stepping by  
     203      !!              applying Asselin time filter and swapping the arrays 
     204      !!              ssha  already computed in ssh_wzv   
     205      !! 
     206      !! ** Method  : - apply Asselin time fiter to now ssh and swap : 
     207      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 
     208      !!             sshn = ssha 
     209      !! 
     210      !! ** action  : - sshb, sshn   : before & now sea surface height 
     211      !!                               ready for the next time step 
     212      !!---------------------------------------------------------------------- 
     213      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     214      !! 
     215      INTEGER  ::   ji, jj               ! dummy loop indices 
     216      !!---------------------------------------------------------------------- 
     217 
     218      IF( kt == nit000 ) THEN 
     219         IF(lwp) WRITE(numout,*) 
     220         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' 
     221         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     222      ENDIF 
     223 
     224      ! Time filter and swap of the ssh 
     225      ! ------------------------------- 
     226 
     227      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points 
     228         !                   ! ---------------------- ! 
     229         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter 
     230            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now) 
     231            sshu_n(:,:) = sshu_a(:,:) 
     232            sshv_n(:,:) = sshv_a(:,:) 
     233            sshf_n(:,:) = sshf_a(:,:) 
     234         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap 
     235            DO jj = 1, jpj 
     236               DO ji = 1, jpi                                ! before <-- now filtered 
     237                  sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) ) 
     238                  sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 
     239                  sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 
     240                  sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 
     241                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after 
     242                  sshu_n(ji,jj) = sshu_a(ji,jj) 
     243                  sshv_n(ji,jj) = sshv_a(ji,jj) 
     244                  sshf_n(ji,jj) = sshf_a(ji,jj) 
     245               END DO 
     246            END DO 
     247         ENDIF 
     248         ! 
     249      ELSE                   ! fixed levels :   ssh at t-point only 
     250         !                   ! ------------ ! 
     251         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter 
     252            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now) 
     253            ! 
     254         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap 
     255            DO jj = 1, jpj 
     256               DO ji = 1, jpi                                ! before <-- now filtered 
     257                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )     
     258                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after 
     259               END DO 
     260            END DO 
     261         ENDIF 
     262         ! 
     263      ENDIF 
     264 
     265      !                      ! write filtered free surface arrays in restart file 
     266      IF( lrst_oce )   CALL ssh_rst( kt, 'WRITE' ) 
     267      ! 
     268      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
     269      ! 
     270   END SUBROUTINE ssh_nxt 
     271 
     272 
     273   SUBROUTINE ssh_rst( kt, cdrw ) 
     274      !!--------------------------------------------------------------------- 
     275      !!                   ***  ROUTINE ssh_rst  *** 
     276      !! 
     277      !! ** Purpose :   Read or write Sea Surface Height arrays in restart file 
     278      !! 
     279      !! ** action  : - cdrw = READ  :  sshb, sshn  read    in ocean restart file 
     280      !!              - cdrw = WRITE :  sshb, sshn  written in ocean restart file 
     281      !!---------------------------------------------------------------------- 
     282      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
     283      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
     284      !!---------------------------------------------------------------------- 
     285      ! 
     286      IF( TRIM(cdrw) == 'READ' ) THEN 
     287         IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
     288            CALL iom_get( numror, jpdom_autoglo, 'sshb'  , sshb(:,:)   ) 
     289            CALL iom_get( numror, jpdom_autoglo, 'sshn'  , sshn(:,:)   ) 
     290            IF( neuler == 0 )   sshb(:,:) = sshn(:,:) 
     291         ELSE 
     292            IF( nn_rstssh == 1 ) THEN 
     293               sshb(:,:) = 0.e0 
     294               sshn(:,:) = 0.e0 
     295            ENDIF 
     296         ENDIF 
     297      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
     298         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb(:,:) ) 
     299         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn(:,:) ) 
     300      ENDIF 
     301      ! 
     302   END SUBROUTINE ssh_rst 
    152303 
    153304   !!====================================================================== 

  • trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90

    r1146 r1438  
    9494            ! of the square root of the resulting N^2 ( required to compute 
    9595            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
    96             zn2 = MAX( rn2(ji,1,jk), 0.e0 ) 
     96            zn2 = MAX( rn2b(ji,1,jk), 0.e0 ) 
    9797            zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 
    9898            ! Compute elements required for the inverse time scale of baroclinic 
    … …  
    113113               ! of the square root of the resulting N^2 ( required to compute  
    114114               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
    115                zn2 = MAX( rn2(ji,jj,jk), 0.e0 ) 
     115               zn2 = MAX( rn2b(ji,jj,jk), 0.e0 ) 
    116116               zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 
    117117               ! Compute elements required for the inverse time scale of baroclinic 

  • trunk/NEMO/OPA_SRC/TRA/tradmp.F90

    r1415 r1438  
    331331      INTEGER ::   ji, jj, jk                   ! dummy loop indices 
    332332      INTEGER ::   ii0, ii1, ij0, ij1           !    "          " 
    333       INTEGER ::   idmp                         ! logical unit for file restoring damping term 
     333      INTEGER ::   inum0                        ! logical unit for file restoring damping term 
    334334      INTEGER ::   icot                         ! logical unit for file distance to the coast 
    335335      REAL(wp) ::   zinfl, zlon                 ! temporary scalars 

  • trunk/NEMO/OPA_SRC/TRA/tranxt.F90

    r1146 r1438  
    1414   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation 
    1515   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf 
     16   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option 
    1617   !!---------------------------------------------------------------------- 
    1718 
    1819   !!---------------------------------------------------------------------- 
    1920   !!   tra_nxt       : time stepping on temperature and salinity 
     21   !!   tra_nxt_fix   : time stepping on temperature and salinity : fixed    volume case 
     22   !!   tra_nxt_vvl   : time stepping on temperature and salinity : variable volume case 
    2023   !!---------------------------------------------------------------------- 
    2124   USE oce             ! ocean dynamics and tracers variables 
    2225   USE dom_oce         ! ocean space and time domain variables  
    2326   USE zdf_oce         ! ??? 
     27   USE domvvl          ! variable volume 
    2428   USE dynspg_oce      ! surface pressure gradient variables 
    2529   USE trdmod_oce      ! ocean variables trends 
    … …  
    3943   PUBLIC   tra_nxt    ! routine called by step.F90 
    4044 
     45   REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 
     46 
    4147   !! * Substitutions 
    4248#  include "domzgr_substitute.h90" 
    … …  
    6773      !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
    6874      !! 
    69       !!              - Apply the Asselin time filter on now fields, 
    70       !!             save in (ta,sa) an average over the three time levels  
    71       !!             which will be used to compute rdn and thus the semi-implicit 
    72       !!             hydrostatic pressure gradient (ln_dynhpg_imp = T), and 
    73       !!             swap tracer fields to prepare the next time_step. 
    74       !!                This can be summurized for tempearture as: 
    75       !!                    zt = (ta+2tn+tb)/4       ln_dynhpg_imp = T 
    76       !!                    zt = 0                   otherwise 
    77       !!                    tb = tn + atfp*[ tb - 2 tn + ta ] 
    78       !!                    tn = ta  
    79       !!                    ta = zt        (NB: reset to 0 after eos_bn2 call) 
     75      !!              - Update lateral boundary conditions on AGRIF children 
     76      !!             domains (lk_agrif=T) 
    8077      !! 
    8178      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
    … …  
    8784      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    8885      !! 
    89       INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    90       REAL(wp) ::   zt, zs, zfact ! temporary scalars 
     86      INTEGER  ::   jk    ! dummy loop indices 
     87      REAL(wp) ::   zfact ! temporary scalars 
    9188      !!---------------------------------------------------------------------- 
    9289 
    … …  
    111108      CALL Agrif_tra                   ! AGRIF zoom boundaries 
    112109#endif 
     110  
     111      ! set time step size (Euler/Leapfrog) 
     112      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler) 
     113      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog) 
     114      ENDIF 
    113115 
    114116      ! trends computation initialisation 
    … …  
    118120      ENDIF 
    119121 
    120       ! Asselin time filter and swap of arrays 
    121       !                                              
    122       IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler 1st time step : swap only 
    123          DO jk = 1, jpkm1 
    124             tb(:,:,jk) = tn(:,:,jk)                       ! ta, sa remain at their values which 
    125             sb(:,:,jk) = sn(:,:,jk)                       ! correspond to tn, sn after the sawp 
    126             tn(:,:,jk) = ta(:,:,jk)                      
    127             sn(:,:,jk) = sa(:,:,jk) 
    128          END DO 
    129          !                                            
    130       ELSE                                           ! Leap-Frog : filter + swap 
    131          !                                             
    132          IF( ln_dynhpg_imp ) THEN                         ! semi-implicite hpg case 
    133             DO jk = 1, jpkm1                              ! (save the averaged of the 3 time steps  
    134                DO jj = 1, jpj                             !                   in the after fields) 
    135                   DO ji = 1, jpi 
    136                      zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25 
    137                      zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25 
    138                      tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 
    139                      sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 
    140                      tn(ji,jj,jk) = ta(ji,jj,jk) 
    141                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    142                      ta(ji,jj,jk) = zt 
    143                      sa(ji,jj,jk) = zs 
    144                   END DO 
    145                END DO 
    146             END DO 
    147          ELSE                                            ! explicit hpg case 
    148             DO jk = 1, jpkm1 
    149                DO jj = 1, jpj 
    150                   DO ji = 1, jpi 
    151                      tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 
    152                      sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 
    153                      tn(ji,jj,jk) = ta(ji,jj,jk) 
    154                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    155                   END DO 
    156                END DO 
    157             END DO 
    158          ENDIF 
    159          ! 
     122      ! Leap-Frog + Asselin filter time stepping 
     123      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt )      ! variable volume level (vvl) 
     124      ELSE                  ;   CALL tra_nxt_fix( kt )      ! fixed    volume level 
    160125      ENDIF 
    161126 
    … …  
    165130#endif       
    166131 
    167       ! trends diagnostics :  Asselin filter trend : (tb filtered - tb)/2dt 
    168       IF( l_trdtra ) THEN      
     132      ! trends computation 
     133      IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt      
    169134         DO jk = 1, jpkm1 
    170             zfact = 1.e0 / ( 2.*rdttra(jk) )             ! NB: euler case, (tb filtered - tb)=0 so 2dt always OK 
     135            zfact = 1.e0 / r2dt_t(jk)              
    171136            ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact 
    172137            ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact 
    … …  
    174139         CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt ) 
    175140      END IF 
     141 
    176142      !                        ! control print 
    177143      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   & 
    … …  
    180146   END SUBROUTINE tra_nxt 
    181147 
     148 
     149   SUBROUTINE tra_nxt_fix( kt ) 
     150      !!---------------------------------------------------------------------- 
     151      !!                   ***  ROUTINE tra_nxt_fix  *** 
     152      !! 
     153      !! ** Purpose :   fixed volume: apply the Asselin time filter and  
     154      !!                swap the tracer fields. 
     155      !!  
     156      !! ** Method  : - Apply a Asselin time filter on now fields. 
     157      !!              - save in (ta,sa) an average over the three time levels  
     158      !!             which will be used to compute rdn and thus the semi-implicit 
     159      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T) 
     160      !!              - swap tracer fields to prepare the next time_step. 
     161      !!                This can be summurized for tempearture as: 
     162      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T 
     163      !!             ztm = 0                   otherwise 
     164      !!             tb  = tn + atfp*[ tb - 2 tn + ta ] 
     165      !!             tn  = ta  
     166      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call) 
     167      !! 
     168      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
     169      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     170      !!---------------------------------------------------------------------- 
     171      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     172      !! 
     173      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     174      REAL(wp) ::   ztm, ztf     ! temporary scalars 
     175      REAL(wp) ::   zsm, zsf     !    -         - 
     176      !!---------------------------------------------------------------------- 
     177 
     178      IF( kt == nit000 ) THEN 
     179         IF(lwp) WRITE(numout,*) 
     180         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping' 
     181         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     182      ENDIF 
     183      ! 
     184      !                                              ! ----------------------- ! 
     185      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case ! 
     186         !                                           ! ----------------------- ! 
     187         ! 
     188         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
     189            DO jk = 1, jpkm1                              ! (only swap) 
     190               DO jj = 1, jpj 
     191                  DO ji = 1, jpi 
     192                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tb <-- tn 
     193                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     194                  END DO 
     195               END DO 
     196            END DO 
     197         ELSE                                             ! general case (Leapfrog + Asselin filter 
     198            DO jk = 1, jpkm1 
     199               DO jj = 1, jpj 
     200                  DO ji = 1, jpi 
     201                     ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! mean t 
     202                     zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 
     203                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t  
     204                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 
     205                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn  
     206                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 
     207                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta 
     208                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     209                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t 
     210                     sa(ji,jj,jk) = zsm 
     211                  END DO 
     212               END DO 
     213            END DO 
     214         ENDIF 
     215         !                                           ! ----------------------- ! 
     216      ELSE                                           !    explicit hpg case    ! 
     217         !                                           ! ----------------------- ! 
     218         ! 
     219         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
     220            DO jk = 1, jpkm1 
     221               DO jj = 1, jpj 
     222                  DO ji = 1, jpi 
     223                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta 
     224                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     225                  END DO 
     226               END DO 
     227            END DO 
     228         ELSE                                             ! general case (Leapfrog + Asselin filter 
     229            DO jk = 1, jpkm1 
     230               DO jj = 1, jpj 
     231                  DO ji = 1, jpi 
     232                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t  
     233                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 
     234                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn  
     235                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 
     236                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta 
     237                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     238                  END DO 
     239               END DO 
     240            END DO 
     241         ENDIF 
     242      ENDIF 
     243      ! 
     244   END SUBROUTINE tra_nxt_fix 
     245 
     246 
     247   SUBROUTINE tra_nxt_vvl( kt ) 
     248      !!---------------------------------------------------------------------- 
     249      !!                   ***  ROUTINE tra_nxt_vvl  *** 
     250      !! 
     251      !! ** Purpose :   Time varying volume: apply the Asselin time filter   
     252      !!                and swap the tracer fields. 
     253      !!  
     254      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
     255      !!              - save in (ta,sa) a thickness weighted average over the three  
     256      !!             time levels which will be used to compute rdn and thus the semi- 
     257      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 
     258      !!              - swap tracer fields to prepare the next time_step. 
     259      !!                This can be summurized for tempearture as: 
     260      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T 
     261      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )     
     262      !!             ztm = 0                                otherwise 
     263      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
     264      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] ) 
     265      !!             tn  = ta  
     266      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call) 
     267      !! 
     268      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
     269      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     270      !!---------------------------------------------------------------------- 
     271      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     272      !!      
     273      INTEGER  ::   ji, jj, jk             ! dummy loop indices 
     274      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar 
     275      REAL(wp) ::   zsm , zsc_f , zsf , zsca, zscn, zscb   !    -         - 
     276      REAL(wp) ::   ze3mr, ze3fr                           !    -         - 
     277      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         - 
     278      !!---------------------------------------------------------------------- 
     279 
     280      IF( kt == nit000 ) THEN 
     281         IF(lwp) WRITE(numout,*) 
     282         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping' 
     283         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     284      ENDIF 
     285 
     286      !                                              ! ----------------------- ! 
     287      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case ! 
     288         !                                           ! ----------------------- ! 
     289         ! 
     290         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
     291            DO jk = 1, jpkm1                              ! (only swap) 
     292               DO jj = 1, jpj 
     293                  DO ji = 1, jpi 
     294                     tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta 
     295                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     296                  END DO 
     297               END DO 
     298            END DO 
     299         ELSE 
     300            DO jk = 1, jpkm1 
     301               DO jj = 1, jpj 
     302                  DO ji = 1, jpi 
     303                     ze3t_b = fse3t_b(ji,jj,jk) 
     304                     ze3t_n = fse3t_n(ji,jj,jk) 
     305                     ze3t_a = fse3t_a(ji,jj,jk) 
     306                     !                                         ! tracer content at Before, now and after 
     307                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b 
     308                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n 
     309                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a 
     310                     ! 
     311                     !                                         ! Asselin filter on thickness and tracer content 
     312                     ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 
     313                     ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   )  
     314                     zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   )  
     315                     ! 
     316                     !                                         ! filtered tracer including the correction  
     317                     ze3fr = 1.e0  / ( ze3t_n + ze3t_f ) 
     318                     ztf   = ze3fr * ( ztcn   + ztc_f  ) 
     319                     zsf   = ze3fr * ( zscn   + zsc_f  ) 
     320                     !                                         ! mean thickness and tracer 
     321                     ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b ) 
     322                     ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   ) 
     323                     zsm   = ze3mr * ( zsca   + 2.* zscn   + zscb   ) 
     324!!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !! 
     325!!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr 
     326                     !                                         ! swap of arrays 
     327                     tb(ji,jj,jk) = ztf                             ! tb <-- tn + filter 
     328                     sb(ji,jj,jk) = zsf 
     329                     tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta 
     330                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     331                     ta(ji,jj,jk) = ztm                             ! ta <-- mean t 
     332                     sa(ji,jj,jk) = zsm 
     333                  END DO 
     334               END DO 
     335            END DO 
     336         ENDIF 
     337         !                                           ! ----------------------- ! 
     338      ELSE                                           !    explicit hpg case    ! 
     339         !                                           ! ----------------------- ! 
     340         ! 
     341         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step 
     342            DO jk = 1, jpkm1                              ! No filter nor thickness weighting computation required 
     343               DO jj = 1, jpj                             ! ONLY swap 
     344                  DO ji = 1, jpi 
     345                     tn(ji,jj,jk) = ta(ji,jj,jk)                                 ! tn <-- ta 
     346                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     347                  END DO 
     348               END DO 
     349            END DO 
     350            !                                             ! general case (Leapfrog + Asselin filter) 
     351         ELSE                                             ! apply filter on thickness weighted tracer and swap 
     352            DO jk = 1, jpkm1 
     353               DO jj = 1, jpj 
     354                  DO ji = 1, jpi 
     355                     ze3t_b = fse3t_b(ji,jj,jk) 
     356                     ze3t_n = fse3t_n(ji,jj,jk) 
     357                     ze3t_a = fse3t_a(ji,jj,jk) 
     358                     !                                         ! tracer content at Before, now and after 
     359                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b 
     360                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n 
     361                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a 
     362                     ! 
     363                     !                                         ! Asselin filter on thickness and tracer content 
     364                     ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 
     365                     ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   )  
     366                     zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   )  
     367                     ! 
     368                     !                                         ! filtered tracer including the correction  
     369                     ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 
     370                     ztf   =  ( ztcn  + ztc_f ) * ze3fr 
     371                     zsf   =  ( zscn  + zsc_f ) * ze3fr 
     372                     !                                         ! swap of arrays 
     373                     tb(ji,jj,jk) = ztf                             ! tb <-- tn filtered 
     374                     sb(ji,jj,jk) = zsf 
     375                     tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta 
     376                     sn(ji,jj,jk) = sa(ji,jj,jk) 
     377                  END DO 
     378               END DO 
     379            END DO 
     380         ENDIF 
     381      ENDIF 
     382      ! 
     383   END SUBROUTINE tra_nxt_vvl 
     384 
    182385   !!====================================================================== 
    183386END MODULE tranxt 

  • trunk/NEMO/OPA_SRC/TRA/trazdf.F90

    r1239 r1438  
    7878         ztrds(:,:,:) = sa(:,:,:) 
    7979      ENDIF 
    80  
    81       IF( lk_vvl )   CALL dom_vvl_ssh( kt )      ! compute ssha field 
    8280 
    8381      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend 

  • trunk/NEMO/OPA_SRC/TRA/trazdf_exp.F90

    r1146 r1438  
    4242#  include "vectopt_loop_substitute.h90" 
    4343   !!---------------------------------------------------------------------- 
    44    !! NEMO/OPA  3.0 , LOCEAN-IPSL (2008)  
     44   !! NEMO/OPA  3.2 , LOCEAN-IPSL (2009)  
    4545   !! $Id$  
    4646   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    7878      INTEGER  ::   ji, jj, jk, jl            ! dummy loop indices 
    7979      REAL(wp) ::   zlavmr, zave3r, ze3tr     ! temporary scalars 
    80       REAL(wp) ::   zta, zsa, ze3tb, zcoef    ! temporary scalars 
     80      REAL(wp) ::   zta, zsa, ze3tb           ! temporary scalars 
    8181      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz, zww   ! 3D workspace 
    8282      !!--------------------------------------------------------------------- 
    … …  
    105105            DO jj = 2, jpjm1  
    106106               DO ji = fs_2, fs_jpim1   ! vector opt. 
    107                   zave3r = 1.e0 / fse3w(ji,jj,jk)  
     107                  zave3r = 1.e0 / fse3w_n(ji,jj,jk)  
    108108                  zwy(ji,jj,jk) =   avt(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) * zave3r 
    109109                  zww(ji,jj,jk) = fsavs(ji,jj,jk) * ( zwz(ji,jj,jk-1) - zwz(ji,jj,jk) ) * zave3r 
    … …  
    115115            DO jj = 2, jpjm1  
    116116               DO ji = fs_2, fs_jpim1   ! vector opt. 
    117                   ze3tr = zlavmr / fse3t(ji,jj,jk) 
     117                  ze3tr = zlavmr / fse3t_n(ji,jj,jk) 
    118118                  zwx(ji,jj,jk) = zwx(ji,jj,jk) + p2dt(jk) * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) * ze3tr 
    119119                  zwz(ji,jj,jk) = zwz(ji,jj,jk) + p2dt(jk) * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) * ze3tr 
    … …  
    130130            DO jj = 2, jpjm1  
    131131               DO ji = fs_2, fs_jpim1   ! vector opt. 
    132                   ze3tb = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )       ! before e3t 
     132                  ze3tb = fse3t_b(ji,jj,jk) / fse3t(ji,jj,jk)                          ! before e3t 
    133133                  zta   = zwx(ji,jj,jk) - tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk)       ! total trends * 2*rdt 
    134134                  zsa   = zwz(ji,jj,jk) - sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk)      
    135                   zcoef = 1.e0 / fse3t(ji,jj,jk) * tmask(ji,jj,jk) 
    136                   ta(ji,jj,jk) = (  ze3tb * tb(ji,jj,jk) + fse3t(ji,jj,jk) * zta  ) * zcoef 
    137                   sa(ji,jj,jk) = (  ze3tb * sb(ji,jj,jk) + fse3t(ji,jj,jk) * zsa  ) * zcoef 
     135                  ta(ji,jj,jk) = (  ze3tb * tb(ji,jj,jk) + zta  ) * tmask(ji,jj,jk) 
     136                  sa(ji,jj,jk) = (  ze3tb * sb(ji,jj,jk) + zsa  ) * tmask(ji,jj,jk) 
    138137               END DO 
    139138            END DO 
    … …  
    143142            DO jj = 2, jpjm1  
    144143               DO ji = fs_2, fs_jpim1   ! vector opt. 
    145                   ta(ji,jj,jk) = ( zwx(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ) *tmask(ji,jj,jk) 
    146                   sa(ji,jj,jk) = ( zwz(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ) *tmask(ji,jj,jk) 
     144                  ta(ji,jj,jk) = ( zwx(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ) * tmask(ji,jj,jk) 
     145                  sa(ji,jj,jk) = ( zwz(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ) * tmask(ji,jj,jk) 
    147146               END DO 
    148147            END DO 

  • trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r1156 r1438  
    11MODULE trazdf_imp 
    2    !!============================================================================== 
     2   !!====================================================================== 
    33   !!                 ***  MODULE  trazdf_imp  *** 
    44   !! Ocean active tracers:  vertical component of the tracer mixing trend 
    5    !!============================================================================== 
    6    !! History : 
    7    !!   6.0  !  90-10  (B. Blanke)  Original code 
    8    !!   7.0  !  91-11 (G. Madec) 
    9    !!        !  92-06 (M. Imbard) correction on tracer trend loops 
    10    !!        !  96-01 (G. Madec) statement function for e3 
    11    !!        !  97-05 (G. Madec) vertical component of isopycnal 
    12    !!        !  97-07 (G. Madec) geopotential diffusion in s-coord 
    13    !!        !  00-08 (G. Madec) double diffusive mixing 
    14    !!   8.5  !  02-08 (G. Madec)  F90: Free form and module 
    15    !!   9.0  !  06-11 (G. Madec)  New step reorganisation 
     5   !!====================================================================== 
     6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code 
     7   !!            7.0  !  1991-11  (G. Madec) 
     8   !!                 !  1992-06  (M. Imbard) correction on tracer trend loops 
     9   !!                 !  1996-01  (G. Madec) statement function for e3 
     10   !!                 !  1997-05  (G. Madec) vertical component of isopycnal 
     11   !!                 !  1997-07  (G. Madec) geopotential diffusion in s-coord 
     12   !!                 !  2000-08  (G. Madec) double diffusive mixing 
     13   !!   NEMO     1.0  !  2002-08  (G. Madec) F90: Free form and module 
     14   !!            2.0  !  2006-11  (G. Madec) New step reorganisation 
     15   !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends 
     16   !!---------------------------------------------------------------------- 
     17   
    1618   !!---------------------------------------------------------------------- 
    1719   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical   
    1820   !!                 part of the mixing tensor. 
    1921   !!---------------------------------------------------------------------- 
    20    !! * Modules used 
    2122   USE oce             ! ocean dynamics and tracers variables 
    2223   USE dom_oce         ! ocean space and time domain variables  
    … …  
    3637   PRIVATE 
    3738 
    38    !! * Routine accessibility 
    39    PUBLIC tra_zdf_imp   !  routine called by step.F90 
     39   PUBLIC   tra_zdf_imp   !  routine called by step.F90 
    4040 
    4141   !! * Substitutions 
    … …  
    4545#  include "vectopt_loop_substitute.h90" 
    4646   !!---------------------------------------------------------------------- 
    47    !!---------------------------------------------------------------------- 
    48    !!  OPA 9.0 , LOCEAN-IPSL (2005)  
     47   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    4948   !! $Id$ 
    5049   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    9089      !! 
    9190      !!--------------------------------------------------------------------- 
    92       !! * Modules used 
    93       USE oce    , ONLY :   zwd   => ua,  &  ! ua used as workspace 
    94                             zws   => va      ! va   "          " 
    95       !! * Arguments 
    96       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index 
    97       REAL(wp), DIMENSION(jpk), INTENT( in ) ::   & 
    98          p2dt                                ! vertical profile of tracer time-step 
    99  
    100       !! * Local declarations 
    101       INTEGER  ::   ji, jj, jk               ! dummy loop indices 
    102       REAL(wp) ::   zavi, zrhs, znvvl,     & ! temporary scalars 
    103          ze3tb, ze3tn, ze3ta, zvsfvvl        ! variable vertical scale factors 
    104       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    105          zwi, zwt, zavsi                     ! workspace arrays 
     91      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace 
     92      USE oce    , ONLY :   zws   => va   ! va  -          - 
     93      !! 
     94      INTEGER                 , INTENT(in) ::   kt     ! ocean time-step index 
     95      REAL(wp), DIMENSION(jpk), INTENT(in) ::   p2dt   ! vertical profile of tracer time-step 
     96      !! 
     97      INTEGER  ::   ji, jj, jk            ! dummy loop indices 
     98      REAL(wp) ::   zavi, zrhs, znvvl     ! temporary scalars 
     99      REAL(wp) ::   ze3tb, ze3tn, ze3ta   ! variable vertical scale factors 
     100      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt, zavsi   ! workspace arrays 
    106101      !!--------------------------------------------------------------------- 
    107102 
    … …  
    169164         DO jj = 2, jpjm1 
    170165            DO ji = fs_2, fs_jpim1   ! vector opt. 
    171                zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) ) 
    172                ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                ! after scale factor at T-point 
    173                ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl                        ! now   scale factor at T-point 
     166               ze3ta =  ( 1. - znvvl )   &                            ! after scale factor at T-point 
     167                  &   +        znvvl    * fse3t_a(ji,jj,jk)  
     168               ze3tn =    znvvl          &                            ! now   scale factor at T-point 
     169                  &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 
    174170               zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) ) 
    175171               zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 
    … …  
    182178      DO jj = 2, jpjm1 
    183179         DO ji = fs_2, fs_jpim1   ! vector opt. 
    184             zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) ) 
    185             ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                   ! after scale factor at T-point 
     180            ze3ta = ( 1. - znvvl )        &                           ! after scale factor at T-point 
     181               &   +       znvvl      * fse3t_a(ji,jj,1)  
    186182            zwi(ji,jj,1) = 0.e0 
    187183            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 
    … …  
    227223      DO jj = 2, jpjm1 
    228224         DO ji = fs_2, fs_jpim1 
    229             zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) ) 
    230             ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl 
    231             ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,1) 
     225            ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1) 
     226            ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1) 
    232227            ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1) 
    233228         END DO 
    … …  
    236231         DO jj = 2, jpjm1 
    237232            DO ji = fs_2, fs_jpim1 
    238                zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) ) 
    239                ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl 
    240                ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,jk) 
     233               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk) 
     234               ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk) 
    241235               zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk)   ! zrhs=right hand side  
    242                ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1) 
     236               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1) 
    243237            END DO 
    244238         END DO 
    … …  
    271265         DO jj = 2, jpjm1 
    272266            DO ji = fs_2, fs_jpim1   ! vector opt. 
    273                zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) ) 
    274                ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                      ! after scale factor at T-point 
    275                ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl                              ! now   scale factor at T-point 
     267               ze3ta =  ( 1. - znvvl )                        &         ! after scale factor at T-point 
     268                  &   +        znvvl   * fse3t_a(ji,jj,jk)            
     269               ze3tn =    znvvl                               &         ! now   scale factor at T-point 
     270                  &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 
    276271               zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) ) 
    277272               zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 
    … …  
    284279      DO jj = 2, jpjm1 
    285280         DO ji = fs_2, fs_jpim1   ! vector opt. 
    286             zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) ) 
    287             ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                          ! after scale factor at T-point 
     281            ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) 
    288282            zwi(ji,jj,1) = 0.e0 
    289283            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 
    … …  
    328322      DO jj = 2, jpjm1 
    329323         DO ji = fs_2, fs_jpim1 
    330             zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) ) 
    331             ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl                               ! before scale factor at T-point 
    332             ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,1)                        ! now    scale factor at T-point 
     324            ze3tb = ( 1. - znvvl )   &                                           ! before scale factor at T-point 
     325               &   +  znvvl       * fse3t_b(ji,jj,1) 
     326            ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,1)                    ! now    scale factor at T-point 
    333327            sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1) 
    334328         END DO 
    … …  
    337331         DO jj = 2, jpjm1 
    338332            DO ji = fs_2, fs_jpim1 
    339                zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) ) 
    340                ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl                            ! before scale factor at T-point 
    341                ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,jk)                    ! now    scale factor at T-point 
     333               ze3tb = ( 1. - znvvl )   &                                        ! before scale factor at T-point 
     334                  &   +  znvvl       * fse3t_b(ji,jj,jk) 
     335               ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)                ! now    scale factor at T-point 
    342336               zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk)     ! zrhs=right hand side 
    343337               sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1) 
    … …  
    361355         END DO 
    362356      END DO 
    363  
     357      ! 
    364358   END SUBROUTINE tra_zdf_imp 
    365359 

  • trunk/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r1152 r1438  
    55   !!                of vertical eddy mixing coefficient 
    66   !!====================================================================== 
     7   !! History :  OPA  !  1997-06  (G. Madec, A. Lazar)  Original code 
     8   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
     9   !!             -   !  2005-06  (C. Ethe) KPP parameterization 
     10   !!            3.2  !  2009-03  (M. Leclair, G. Madec, R. Benshila) test on both before & after 
     11   !!---------------------------------------------------------------------- 
    712 
    813   !!---------------------------------------------------------------------- 
    9    !!   zdf_evd      : update momentum and tracer Kz at the location of 
    10    !!                  statically unstable portion of the water column 
    11    !!                  (called if ln_zdfevd=T) 
     14   !!   zdf_evd      : increase the momentum and tracer Kz at the location of 
     15   !!                  statically unstable portion of the water column (ln_zdfevd=T) 
    1216   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    1417   USE oce             ! ocean dynamics and tracers variables 
    1518   USE dom_oce         ! ocean space and time domain variables 
    … …  
    2225   PRIVATE 
    2326 
    24    !! * Routine accessibility 
    25    PUBLIC zdf_evd      ! called by step.F90 
     27   PUBLIC   zdf_evd    ! called by step.F90 
    2628 
    2729   !! * Substitutions 
    2830#  include "domzgr_substitute.h90" 
    2931   !!---------------------------------------------------------------------- 
    30    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
    31    !! $Id$  
    32    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     32   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
     33   !! $Id$ 
     34   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
    3335   !!---------------------------------------------------------------------- 
    3436 
    … …  
    4850      !! ** Action  :   Update avt, avmu, avmv in statically instable cases 
    4951      !!                and avt_evd which is avt due to convection 
    50       !! References : 
    51       !!      Lazar, A., these de l'universite Paris VI, France, 1997 
    52       !! History : 
    53       !!   7.0  !  97-06  (G. Madec, A. Lazar)  Original code 
    54       !!   8.5  !  02-06  (G. Madec)  F90: Free form and module 
    55       !!   9.0  !  05-06  (C. Ethe) KPP parameterization 
     52      !! References :   Lazar, A., these de l'universite Paris VI, France, 1997 
    5653      !!---------------------------------------------------------------------- 
    57       !! * Arguments 
    5854      INTEGER, INTENT( in ) ::   kt         ! ocean time-step indexocean time step 
    59  
    60       !! * Local declarations 
     55      !! 
    6156      INTEGER ::   ji, jj, jk               ! dummy loop indices 
    6257      !!---------------------------------------------------------------------- 
    … …  
    7974         DO jk = 1, jpkm1                                 ! Horizontal slab 
    8075            !                                             ! =============== 
    81 #   if defined key_vectopt_loop 
    82 !!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
    83             jj = 1                     ! big loop forced 
    84             DO ji = jpi+2, jpij    
    85 #   if defined key_zdfkpp 
    86 !! no implicit mixing in the boundary layer with KPP 
    87                IF( ( rn2(ji,jj,jk) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN 
    88 #   else 
    89                IF( rn2(ji,jj,jk) <= -1.e-12 ) THEN 
    90 #   endif 
    91                   avt (ji  ,jj  ,jk) = avevd * tmask(ji  ,jj  ,jk) 
    92                   avmu(ji  ,jj  ,jk) = avevd * umask(ji  ,jj  ,jk) 
    93                   avmu(ji-1,jj  ,jk) = avevd * umask(ji-1,jj  ,jk) 
    94                   avmv(ji  ,jj  ,jk) = avevd * vmask(ji  ,jj  ,jk) 
    95                   avmv(ji  ,jj-1,jk) = avevd * vmask(ji  ,jj-1,jk) 
    96                ENDIF 
    97             END DO 
    98 #   else 
     76#if defined key_vectopt_loop 
     77            DO jj = 1, 1                     ! big loop forced 
     78               DO ji = jpi+2, jpij    
     79#else 
    9980            DO jj = 2, jpj             ! no vector opt. 
    10081               DO ji = 2, jpi 
    101 #   if defined key_zdfkpp 
     82#endif 
     83#if defined key_zdfkpp 
    10284!! no implicit mixing in the boundary layer with KPP 
    103                IF( ( rn2(ji,jj,jk) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN 
    104 #   else 
    105                IF( rn2(ji,jj,jk) <= -1.e-12 ) THEN 
    106 #   endif 
     85                  IF( ( MIN( rn2(ji,jj,jk),  rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN 
     86#else 
     87                  IF(   MIN( rn2(ji,jj,jk),  rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 
     88#endif 
    10789                     avt (ji  ,jj  ,jk) = avevd * tmask(ji  ,jj  ,jk) 
    10890                     avmu(ji  ,jj  ,jk) = avevd * umask(ji  ,jj  ,jk) 
    … …  
    11395               END DO 
    11496            END DO 
    115 #   endif 
    11697            !                                             ! =============== 
    11798         END DO                                           !   End of slab 
    … …  
    129110            !                                             ! =============== 
    130111!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL!  
    131 #   if defined key_vectopt_loop 
    132             jj = 1                     ! big loop forced 
    133             DO ji = 1, jpij    
    134 #   if defined key_zdfkpp 
    135 !! no implicit mixing in the boundary layer with KPP 
    136                IF( ( rn2(ji,jj,jk) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) &               
    137                   avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 
    138 #   else 
    139                IF( rn2(ji,jj,jk) <= -1.e-12 )   avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 
    140 #   endif 
    141             END DO 
    142 #   else 
     112#if defined key_vectopt_loop 
     113            DO jj = 1, 1                     ! big loop forced 
     114               DO ji = 1, jpij    
     115#else 
    143116            DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call) 
    144117               DO ji = 1, jpi 
    145 #   if defined key_zdfkpp 
     118#endif 
     119#if defined key_zdfkpp 
    146120!! no implicit mixing in the boundary layer with KPP 
    147                IF( ( rn2(ji,jj,jk) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) &           
    148                   avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 
    149 #   else 
    150                   IF( rn2(ji,jj,jk) <= -1.e-12 )   avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 
    151 #   endif 
     121                  IF( ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) &           
     122#else 
     123                  IF(   MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) & 
     124#endif 
     125                     avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 
    152126               END DO 
    153127            END DO 
    154 #   endif 
    155128            !                                             ! =============== 
    156129         END DO                                           !   End of slab 
    … …  
    159132 
    160133      ! update of avt_evd and avmu_evd 
    161       avt_evd  (:,:,:) = avt (:,:,:)  - avt_evd (:,:,:)  
    162       avmu_evd (:,:,:) = avmu(:,:,:)  - avmu_evd (:,:,:)  
     134      avt_evd (:,:,:) = avt (:,:,:)  - avt_evd (:,:,:)  
     135      avmu_evd(:,:,:) = avmu(:,:,:)  - avmu_evd(:,:,:)  
    163136 
    164137   END SUBROUTINE zdf_evd 

  • trunk/NEMO/OPA_SRC/ZDF/zdftke.F90

    r1268 r1438  
    224224!CDIR NOVERRCHK 
    225225            DO ji = fs_2, fs_jpim1   ! vector opt. 
    226                zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
     226               zrn2 = MAX( rn2b(ji,jj,jk), rsmall ) 
    227227               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  ) 
    228228            END DO 
    … …  
    321321         ! 
    322322         ! Computation of total energy produce by LC : cumulative sum over jk 
    323          zpelc(:,:,1) =  MAX( rn2(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1) 
     323         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1) 
    324324         DO jk = 2, jpk 
    325             zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 
     325            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 
    326326         END DO 
    327327         ! 
    … …  
    427427                  zdkv = zcoef * (  vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   & 
    428428                     &            - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  ) 
    429                   zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)        ! coefficient (zesh2) 
     429                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2b(ji,jj,jk)        ! coefficient (zesh2) 
    430430                  ! 
    431431                  !                                                             ! Matrix 
    … …  
    464464                  &               - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )  ) 
    465465                  zsh2 = zdku * zdku + zdkv * zdkv                                   ! square of shear 
    466                   zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )                ! local Richardson number 
     466                  zri  = MAX( rn2b(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )                ! local Richardson number 
    467467# if defined key_c1d 
    468468                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri  
    … …  
    471471                  IF( zri >= 0.2 )   zpdl = 0.2 / zri 
    472472                  zpdl = MAX( 0.1, zpdl ) 
    473                   zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)                       ! coefficient (esh2) 
     473                  zesh2 = eboost * zsh2 - zpdl * rn2b(ji,jj,jk)                       ! coefficient (esh2) 
    474474                  ! 
    475475                  !                                                             ! Matrix 

  • trunk/NEMO/OPA_SRC/istate.F90

    r1200 r1438  
    132132          
    133133      ENDIF 
    134  
    135       IF( lk_vvl .OR. lk_agrif ) THEN 
     134      ! 
     135      IF( lk_agrif ) THEN 
    136136         ! read free surface arrays in restart file 
    137137         IF( ln_rstart ) THEN 
    138138            IF( lk_dynspg_flt )   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields 
    139             !                                                         ! gcx, gcxb, sshb, sshn 
    140             IF( lk_dynspg_ts  )   CALL ts_rst ( nit000, 'READ' )      ! read or initialize the following fields 
    141             !                                                         ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 
    142             IF( lk_dynspg_exp )   CALL exp_rst( nit000, 'READ' )      ! read or initialize the following fields 
    143             !                                                         ! sshb, sshn 
    144          ENDIF 
     139            !                                                         ! gcx, gcxb for agrif_opa_init 
     140         ENDIF                                                        ! explicit case not coded yet with AGRIF 
    145141      ENDIF 
    146  
    147       IF( lk_vvl ) THEN 
    148  
    149          ! 
    150          IF( .NOT. lk_dynspg_flt ) sshbb(:,:) = sshb(:,:) 
    151          ! 
    152          CALL dom_vvl               ! ssh init and vertical grid update 
    153  
    154          CALL eos_init              ! usefull to get the equation state type neos parameter 
    155  
    156          CALL bn2( tb, sb, rn2 )    ! before Brunt Vaissala frequency 
    157  
    158          IF( .NOT. ln_rstart ) CALL wzv( nit000 )  
    159  
    160       ENDIF 
    161  
    162       !                                       ! Vertical velocity 
    163       !                                       ! ----------------- 
    164  
    165       IF( .NOT. lk_vvl )    CALL wzv( nit000 )                         ! from horizontal divergence 
    166142      ! 
    167143   END SUBROUTINE istate_init 

  • trunk/NEMO/OPA_SRC/oce.F90

    r1152 r1438  
    44   !! Ocean        :  dynamics and active tracers defined in memory  
    55   !!====================================================================== 
    6    !! History : 
    7    !!   8.5  !  02-11  (G. Madec)  F90: Free form and module 
    8    !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
     6   !! History :  0.1  !  2002-11  (G. Madec)  F90: Free form and module 
     7   !!            1.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
     8   !!            3.1  !  2009-02  (G. Madec, M. Leclair)  pure z* coordinate 
    99   !!---------------------------------------------------------------------- 
    10    !!  OPA 9.0 , LOCEAN-IPSL (2005)  
    11    !! $Id$  
    12    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    13    !!---------------------------------------------------------------------- 
    14    !! * Modules used 
    1510   USE par_oce      ! ocean parameters 
    1611 
    … …  
    2015   !! Physics and algorithm flags 
    2116   !! --------------------------- 
    22    LOGICAL, PUBLIC ::   l_traldf_rot    = .FALSE.  !: rotated laplacian operator for lateral diffusion  
     17   LOGICAL, PUBLIC ::   l_traldf_rot    = .FALSE.  !: rotated laplacian operator for lateral diffusion 
    2318   LOGICAL, PUBLIC ::   ln_dynhpg_imp   = .FALSE.  !: semi-implicite hpg flag 
    2419   INTEGER, PUBLIC ::   nn_dynhpg_rst   = 0        !: add dynhpg implicit variables in restart ot not 
    2520 
    26    !! dynamics and tracer fields 
    27    !! -------------------------- 
    28    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !: 
    29       ! before !  now      !  after  !      ! the after trends becomes the fields 
    30       ! fields !  fields   !  trends !      ! only in dyn(tra)_zdf and dyn(tra)_nxt 
    31       ub       ,  un       ,  ua     ,   &  !: i-horizontal velocity (m/s) 
    32       vb       ,  vn       ,  va     ,   &  !: j-horizontal velocity (m/s) 
    33                   wn       ,             &  !: vertical velocity (m/s) 
    34       rotb     ,  rotn     ,             &  !: relative vorticity (1/s) 
    35       hdivb    ,  hdivn    ,             &  !: horizontal divergence (1/s) 
    36       tb       ,  tn       ,  ta     ,   &  !: potential temperature (celcius) 
    37       sb       ,  sn       ,  sa            !: salinity (psu) 
    38    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !: 
    39       rhd ,                              &  !: in situ density anomalie rhd=(rho-rau0)/rau0 (no units) 
    40       rhop,                              &  !: potential volumic mass (kg/m3) 
    41       rn2                                   !: brunt-vaisala frequency (1/s2) 
     21   !! dynamics and tracer fields             !  before  !  now     !  after   ! the after trends becomes the fields 
     22   !! --------------------------             !  fields  !  fields  !  trends  ! only after tra_zdf and dyn_spg 
     23   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   ub     ,  un      ,  ua      !: i-horizontal velocity      [m/s] 
     24   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vb     ,  vn      ,  va      !: j-horizontal velocity      [m/s] 
     25   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::             wn                 !: vertical velocity          [m/s] 
     26   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rotb   ,  rotn               !: relative vorticity         [s-1] 
     27   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   hdivb  ,  hdivn              !: horizontal divergence      [s-1] 
     28   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tb     ,  tn      ,  ta      !: potential temperature      [Celcius] 
     29   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   sb     ,  sn      ,  sa      !: salinity                   [psu] 
     30   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rn2b   ,  rn2                !: brunt-vaisala frequency**2 [s-2] 
     31   ! 
     32   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0     [no units] 
     33   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rhop   !: potential volumic mass                           [kg/m3] 
    4234 
    43       !! advection scheme choice 
    44       !! ----------------------- 
    45       CHARACTER(len=3), PUBLIC  ::   l_adv   !: 'ce2' centre scheme used 
    46          !                                   !: 'tvd' TVD scheme used 
    47          !                                   !: 'mus' MUSCL scheme used 
    48          !                                   !: 'mu2' MUSCL2 scheme used 
     35   !! advection scheme choice 
     36   !! ----------------------- 
     37   CHARACTER(len=3), PUBLIC  ::   l_adv   !: flag for the advection scheme used (= 'ce2', 'tvd', 'mus' or ...) 
    4938 
    5039   !! surface pressure gradient 
    5140   !! ------------------------- 
    52    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !: 
    53       spgu, spgv             !: horizontal surface pressure gradient 
     41   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   spgu, spgv      !: horizontal surface pressure gradient 
    5442 
    5543   !! interpolated gradient (only used in zps case) 
    5644   !! --------------------- 
    57    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !: 
    58       gtu, gsu, gru,      &  !: t-, s- and rd horizontal gradient at u- and  
    59       gtv, gsv, grv          !: v-points at bottom ocean level  
     45   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtu, gsu, gru   !: horizontal gradient of T, S and rd at bottom u-point 
     46   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtv, gsv, grv   !: horizontal gradient of T, S and rd at bottom v-point  
    6047 
    61    !! free surface 
    62    !! ------------ 
    63    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !: 
    64       sshb , sshn ,        &  !: before, now sea surface height (meters) 
    65       sshu , sshv ,        &  !: sea surface height at u- and v- point 
    66       sshbb, ssha             !: before before sea surface height at t-point 
     48   !! free surface                       !  before  !  now     !  after   ! 
     49   !! ------------                       !  fields  !  fields  !  trends  ! 
     50   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshb   ,  sshn    ,  ssha    !: sea surface height at t-point [m] 
     51   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshu_b ,  sshu_n  ,  sshu_a  !: sea surface height at u-point [m] 
     52   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshv_b ,  sshv_n  ,  sshv_a  !: sea surface height at u-point [m] 
     53   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshf_b ,  sshf_n  ,  sshf_a  !: sea surface height at f-point [m] 
    6754 
    6855#if defined key_dynspg_rl   ||   defined key_esopa 
    6956   !! rigid-lid formulation 
    7057   !! --------------------- 
    71    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !: 
    72       bsfb, bsfn,         &  !: before, now barotropic streamfunction (m3/s) 
    73       bsfd                   !: now trend of barotropic streamfunction (m3/s2) 
     58   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bsfb, bsfn   !: before, now barotropic streamfunction (m3/s) 
     59   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bsfd         !: now trend of barotropic streamfunction (m3/s2) 
    7460#endif 
    7561   !!---------------------------------------------------------------------- 
     62   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2008)  
     63   !! $Id$  
     64   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     65   !!====================================================================== 
    7666END MODULE oce 

  • trunk/NEMO/OPA_SRC/restart.F90

    r1409 r1438  
    143143      CALL iom_rstput( kt, nitrst, numrow, 'rotn'   , rotn    ) 
    144144      CALL iom_rstput( kt, nitrst, numrow, 'hdivn'  , hdivn   ) 
     145      CALL iom_rstput( kt, nitrst, numrow, 'rn2' , rn2  ) 
    145146 
    146147#if defined key_zdftke2 
    … …  
    150151         CALL iom_rstput( kt, nitrst, numrow, 'rhd' , rhd  ) 
    151152#if defined key_zdftke2 
    152          CALL iom_rstput( kt, nitrst, numrow, 'rn2' , rn2  ) 
    153153         CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt  ) 
    154154      ENDIF 
    … …  
    239239         hdivb(:,:,:) = hdivn(:,:,:) 
    240240      ENDIF 
     241      CALL eos_init            ! usefull to get the equation state type neos parameter 
     242      IF( iom_varid( numror, 'rn2', ldstop = .FALSE. ) > 0 ) THEN 
     243         CALL iom_get( numror, jpdom_autoglo, 'rn2' , rn2  ) 
     244      ELSE 
     245         CALL bn2( tb, sb, rn2 )  ! before Brunt-Vaisala frequency    
     246      ENDIF 
    241247 
    242248#if defined key_zdftke2 
    … …  
    257263#endif 
    258264#if defined key_zdftke2 
    259       CALL eos_init            ! usefull to get the equation state type neos parameter 
    260       IF( iom_varid( numror, 'rn2', ldstop = .FALSE. ) > 0 ) THEN 
    261          CALL iom_get( numror, jpdom_autoglo, 'rn2' , rn2  ) 
    262       ELSE 
    263          IF ( ln_dynhpg_imp ) THEN 
    264             CALL bn2( tb, sb, rn2 )  ! before Brunt-Vaisala frequency    
    265          ENDIF       
    266       ENDIF 
    267265      IF( iom_varid( numror, 'avt', ldstop = .FALSE. ) > 0 ) THEN 
    268266         CALL iom_get( numror, jpdom_autoglo, 'avt' , avt  ) 

  • trunk/NEMO/OPA_SRC/step.F90

    r1417 r1438  
    44   !! Time-stepping    : manager of the ocean, tracer and ice time stepping 
    55   !!====================================================================== 
    6    !! History :        !  91-03  (G. Madec)  Original code 
    7    !!                  !  92-06  (M. Imbard)  add a first output record 
    8    !!                  !  96-04  (G. Madec)  introduction of dynspg 
    9    !!                  !  96-04  (M.A. Foujols)  introduction of passive tracer 
    10    !!             8.0  !  97-06  (G. Madec)  new architecture of call 
    11    !!             8.2  !  97-06  (G. Madec, M. Imbard, G. Roullet)  free surface 
    12    !!             8.2  !  99-02  (G. Madec, N. Grima)  hpg implicit 
    13    !!             8.2  !  00-07  (J-M Molines, M. Imbard)  Open Bondary Conditions 
    14    !!             9.0  !  02-06  (G. Madec)  free form, suppress macro-tasking 
    15    !!             " "  !  04-08  (C. Talandier) New trends organization 
    16    !!             " "  !  05-01  (C. Ethe) Add the KPP closure scheme 
    17    !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization 
    18    !!             " "  !  05-11  (G. Madec)  Reorganisation of tra and dyn calls 
    19    !!             " "  !  06-01  (L. Debreu, C. Mazauric)  Agrif implementation 
    20    !!             " "  !  06-07  (S. Masson)  restart using iom 
    21    !!             " "  !  06-08  (G. Madec)  surface module  
    22    !!             " "  !  07-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY) 
     6   !! History :  OPA  !  1991-03  (G. Madec)  Original code 
     7   !!                 !  1991-11  (G. Madec) 
     8   !!                 !  1992-06  (M. Imbard)  add a first output record 
     9   !!                 !  1996-04  (G. Madec)  introduction of dynspg 
     10   !!                 !  1996-04  (M.A. Foujols)  introduction of passive tracer 
     11   !!            8.0  !  1997-06  (G. Madec)  new architecture of call 
     12   !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface 
     13   !!             -   !  1999-02  (G. Madec, N. Grima)  hpg implicit 
     14   !!             -   !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions 
     15   !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking 
     16   !!             -   !  2004-08  (C. Talandier) New trends organization 
     17   !!                 !  2005-01  (C. Ethe) Add the KPP closure scheme 
     18   !!                 !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls 
     19   !!                 !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation 
     20   !!                 !  2006-07  (S. Masson)  restart using iom 
     21   !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate 
    2322   !!---------------------------------------------------------------------- 
    2423 
    … …  
    117116   USE restart         ! ocean restart                    (rst_wri routine) 
    118117   USE prtctl          ! Print control                    (prt_ctl routine) 
    119    USE domvvl          ! variable volume                  (dom_vvl routine) 
    120118 
    121119#if defined key_agrif 
    … …  
    141139#if defined key_agrif 
    142140   SUBROUTINE stp( ) 
     141      INTEGER             ::   kstp   ! ocean time-step index 
    143142#else 
    144143   SUBROUTINE stp( kstp ) 
     144      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
    145145#endif 
    146146      !!---------------------------------------------------------------------- 
    … …  
    160160      !!              -8- Outputs and diagnostics 
    161161      !!---------------------------------------------------------------------- 
    162       !! * Arguments 
    163 #if defined key_agrif    
    164       INTEGER             ::   kstp   ! ocean time-step index 
    165 #else 
    166       INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
    167 #endif       
    168162      INTEGER ::   jk       ! dummy loop indice 
    169163      INTEGER ::   indic    ! error indicator if < 0 
    … …  
    202196      ENDIF 
    203197 
     198      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     199      !                Ocean dynamics : ssh, wn, hdiv, rot                   ! 
     200      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     201                       CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity 
     202      IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence) 
     203                       CALL ssh_wzv( kstp )                 ! after ssh & vertical velocity 
    204204 
    205205      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    206206      ! Ocean physics update 
    207207      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     208      ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
     209      !----------------------------------------------------------------------- 
     210      IF( neuler == 0 .AND. kstp == nit000 ) THEN 
     211                        CALL bn2( tn, sn, rn2 )              ! now    Brunt-Vaisala frequency 
     212                        rn2b(:,:,:) = rn2(:,:,:) 
     213      ELSE 
     214                        rn2b(:,:,:) = rn2(:,:,:)             ! before Brunt-Vaisala frequency   
     215                        CALL bn2( tn, sn, rn2 )              ! now    Brunt-Vaisala frequency 
     216      ENDIF 
     217 
    208218#if defined key_zdftke2 
    209219      IF ( ln_dynhpg_imp ) THEN 
    … …  
    211221      !  LATERAL PHYSICS 
    212222      !----------------------------------------------------------------------- 
    213       ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
    214       !----------------------------------------------------------------------- 
    215223                               CALL zdf_mxl( kstp )                 ! mixed layer depth 
    216          IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2 )       ! before slope of the lateral mixing 
     224         IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing 
    217225#  if defined key_traldf_c2d 
    218226         IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )                 ! eddy induced velocity coefficient 
    … …  
    223231      !  VERTICAL PHYSICS 
    224232      !----------------------------------------------------------------------- 
    225       ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
    226       !----------------------------------------------------------------------- 
    227 #if defined key_zdftke2 
    228                         CALL bn2( tn, sn, rn2 )              ! now Brunt-Vaisala frequency 
    229 #else 
    230                         CALL bn2( tb, sb, rn2 )              ! before Brunt-Vaisala frequency 
    231 #endif 
    232233      !                                                     ! Vertical eddy viscosity and diffusivity coefficients 
    233234      IF( lk_zdfric )   CALL zdf_ric( kstp )                       ! Richardson number dependent Kz 
    … …  
    261262#if defined key_zdftke2 
    262263      IF( .NOT. ln_dynhpg_imp ) THEN 
    263                         CALL bn2( tb, sb, rn2 )              ! before Brunt-Vaisala frequency 
    264264                        CALL eos( tb, sb, rhd, rhop )        ! now (swap=before) in situ density for dynhpg module 
    265265#endif 
    … …  
    267267      !  LATERAL PHYSICS 
    268268      !----------------------------------------------------------------------- 
    269       ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
    270       !----------------------------------------------------------------------- 
    271       IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2 )       ! before slope of the lateral mixing 
     269      IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing 
    272270#if defined key_traldf_c2d 
    273271      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )                 ! eddy induced velocity coefficient 
    … …  
    363361                               CALL dyn_spg( kstp, indic )    ! surface pressure gradient 
    364362                               CALL dyn_nxt( kstp )           ! lateral velocity at next time step 
    365       IF( lk_vvl )             CALL dom_vvl                   ! vertical mesh at next time step 
    366  
    367  
    368       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    369       ! Computation of diagnostic variables 
    370       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    371       ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
    372       !----------------------------------------------------------------------- 
    373                        CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity 
    374       IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence) 
    375                        CALL wzv( kstp )                     ! Vertical velocity 
     363 
     364                               CALL ssh_nxt( kstp )           ! sea surface height at next time step 
    376365 
    377366      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.