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

Changeset 455


Ignore:
Timestamp:
2006-05-10T18:53:54+02:00 (18 years ago)
Author:
opalod
Message:

nemo_v1_update_048:RB: reorganization of dynamics part, in addition change atsk to jki, suppress dynhpg_atsk.F90 dynzdf_imp_atsk.F90 dynzdf_iso.F90

Location:
trunk/NEMO/OPA_SRC/DYN
Files:
3 deleted
17 edited

Legend:

Unmodified
Added
Removed

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

    r392 r455  
    5151      !!      (note that the Asselin filter has not been applied on hdivb) 
    5252      !!         - compute the now divergence given by : 
    53       !!            * s-coordinate ('key_s_coord' defined) 
    5453      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    55       !!         * z-coordinate (default key) 
    56       !!         hdivn = 1/(e1t*e2t) [ di(e2u  un) + dj(e1v  vn) ] 
     54      !!      Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the 
     55      !!      above expression 
    5756      !!         - apply lateral boundary conditions on hdivn  
    5857      !!      II. vorticity : 
    … …  
    109108         DO jj = 2, jpjm1 
    110109            DO ji = fs_2, fs_jpim1   ! vector opt. 
    111 #if defined key_s_coord || defined key_partial_steps 
    112                hdivn(ji,jj,jk) =   & 
    113                   (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk)  * un(ji-1,jj  ,jk)       & 
    114                    + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk)  * vn(ji  ,jj-1,jk)  )    & 
    115                   / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    116 #else 
     110#if defined key_zco 
    117111               hdivn(ji,jj,jk) = (  e2u(ji,jj) * un(ji,jj,jk) - e2u(ji-1,jj  ) * un(ji-1,jj  ,jk)      & 
    118112                  &               + e1v(ji,jj) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * vn(ji  ,jj-1,jk)  )   & 
    119      &            / ( e1t(ji,jj) * e2t(ji,jj) ) 
     113                  &            / ( e1t(ji,jj) * e2t(ji,jj) ) 
     114#else 
     115               hdivn(ji,jj,jk) =   & 
     116                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       & 
     117                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    & 
     118                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    120119#endif 
    121120            END DO 
    … …  
    130129         IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south 
    131130#endif          
    132 #if defined key_agrif 
    133          if ( .NOT. AGRIF_Root() ) then 
    134             IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east 
    135             IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west 
    136             IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north 
    137             IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south 
    138          endif 
    139 #endif        
    140131 
    141132         !                                             ! -------- 
    … …  
    260251      !!      (note that the Asselin filter has not been applied on hdivb) 
    261252      !!      - compute the now divergence given by : 
    262       !!         * s-coordinate ('key_s_coord' defined) 
    263253      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    264       !!         * z-coordinate (default key) 
    265       !!         hdivn = 1/(e1t*e2t) [ di(e2u  un) + dj(e1v  vn) ] 
     254      !!      Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the  
     255      !!      above expression 
    266256      !!      - apply lateral boundary conditions on hdivn  
    267257      !!              - Relavtive Vorticity : 
    … …  
    313303         DO jj = 2, jpjm1 
    314304            DO ji = fs_2, fs_jpim1   ! vector opt. 
    315 #if defined key_s_coord || defined key_partial_steps 
    316                hdivn(ji,jj,jk) =   & 
    317                   (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk)  * un(ji-1,jj  ,jk)       & 
    318                    + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk)  * vn(ji  ,jj-1,jk)  )    & 
    319                   / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    320 #else 
     305#if defined key_zco 
    321306               hdivn(ji,jj,jk) = (  e2u(ji,jj) * un(ji,jj,jk) - e2u(ji-1,jj  ) * un(ji-1,jj  ,jk)      & 
    322307                  &               + e1v(ji,jj) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * vn(ji  ,jj-1,jk)  )   &  
    323308                  / ( e1t(ji,jj) * e2t(ji,jj) ) 
     309#else 
     310               hdivn(ji,jj,jk) =   & 
     311                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       & 
     312                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    & 
     313                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    324314#endif 
    325315            END DO   
    … …  
    334324         IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south 
    335325#endif          
    336 #if defined key_agrif 
    337          if ( .NOT. AGRIF_Root() ) then 
    338             IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east 
    339             IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west 
    340             IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north 
    341             IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south 
    342          endif 
    343 #endif        
    344326         !                                             ! -------- 
    345327         ! relative vorticity                          !   rot  

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

    r359 r455  
    66 
    77   !!---------------------------------------------------------------------- 
    8    !!   dyn_hpg      : update the momentum trend with the horizontal 
     8   !!   dyn_hpg      : update the momentum trend with the now horizontal 
    99   !!                  gradient of the hydrostatic pressure 
    10    !! 
    11    !!   default case : use of 3D work arrays (vector opt. available) 
    12    !!   key_s_coord       : s-coordinate 
    13    !!   key_partial_steps : z-coordinate with partial steps 
    14    !!   default key       : z-coordinate 
     10   !!                  default case : k-j-i loops (vector opt. available) 
     11   !!       hpg_ctl  : initialisation and control of options 
     12   !!       hpg_zco  : z-coordinate scheme 
     13   !!       hpg_zps  : z-coordinate plus partial steps (interpolation) 
     14   !!       hpg_sco  : s-coordinate (standard jacobian formulation) 
     15   !!       hpg_hel  : s-coordinate (helsinki modification) 
     16   !!       hpg_wdj  : s-coordinate (weighted density jacobian) 
     17   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial) 
     18   !!       hpg_rot  : s-coordinate (ROTated axes scheme) 
    1519   !!---------------------------------------------------------------------- 
    1620   !! * Modules used 
    1721   USE oce             ! ocean dynamics and tracers 
    1822   USE dom_oce         ! ocean space and time domain 
     23   USE dynhpg_jki      ! 
    1924   USE phycst          ! physical constants 
    2025   USE in_out_manager  ! I/O manager 
    … …  
    2227   USE trdmod_oce      ! ocean variables trends 
    2328   USE prtctl          ! Print control 
     29   USE lbclnk          ! lateral boundary condition  
    2430 
    2531   IMPLICIT NONE 
    … …  
    2834   !! * Accessibility 
    2935   PUBLIC dyn_hpg                ! routine called by step.F90 
     36 
     37#if defined key_mpp_omp 
     38   !!---------------------------------------------------------------------- 
     39   !!   'key_mpp_omp' :                                 j-k-i loop (j-slab) 
     40   !!---------------------------------------------------------------------- 
     41   LOGICAL, PUBLIC, PARAMETER ::   lk_dynhpg_jki = .TRUE.    !: OpenMP hpg flag 
     42   LOGICAL, PUBLIC, PARAMETER ::   lk_dynhpg     = .FALSE.   !: vector hpg flag 
     43#else 
     44   !!---------------------------------------------------------------------- 
     45   !!   default case :                             k-j-i loop (vector opt.) 
     46   !!----------------------------------------------------------------------    
     47   LOGICAL, PUBLIC, PARAMETER ::   lk_dynhpg_jki = .FALSE.   !: OpenMP hpg flag 
     48   LOGICAL, PUBLIC, PARAMETER ::   lk_dynhpg     = .TRUE.    !: vector hpg flag 
     49#endif 
     50 
     51   !! * Share module variables 
     52   LOGICAL  ::               & !!! ** nam_dynhpg **   hpg flags 
     53      ln_hpg_zco = .TRUE. ,  &  ! z-coordinate - full steps 
     54      ln_hpg_zps = .FALSE.,  &  ! z-coordinate - partial steps (interpolation) 
     55      ln_hpg_sco = .FALSE.,  &  ! s-coordinate (standard jacobian formulation) 
     56      ln_hpg_hel = .FALSE.,  &  ! s-coordinate (helsinki modification) 
     57      ln_hpg_wdj = .FALSE.,  &  ! s-coordinate (weighted density jacobian) 
     58      ln_hpg_djc = .FALSE.,  &  ! s-coordinate (Density Jacobian with Cubic polynomial) 
     59      ln_hpg_rot = .FALSE.      ! s-coordinate (ROTated axes scheme) 
     60 
     61   REAL(wp) ::               & !!! ** nam_dynhpg **  
     62      gamm = 0.e0               ! weighting coefficient 
     63 
     64   INTEGER  ::               &  !  
     65      nhpg  =  0                 ! = 0 to 6, type of pressure gradient scheme used 
     66      !                         ! (deduced from ln_hpg_... flags) 
    3067 
    3168   !! * Substitutions 
    … …  
    4077CONTAINS 
    4178 
    42 #if defined key_s_coord 
    43    !!---------------------------------------------------------------------- 
    44    !!   'key_s_coord' :                                        s-coordinate 
    45    !!----------------------------------------------------------------------    
    46  
    4779   SUBROUTINE dyn_hpg( kt ) 
    4880      !!--------------------------------------------------------------------- 
    4981      !!                  ***  ROUTINE dyn_hpg  *** 
    5082      !! 
    51       !! ** Purpose :   Compute the now momentum trend due to the hor. gradient 
    52       !!      of the hydrostatic pressure. Add it to the general momentum trend. 
    53       !! 
    54       !! ** Method  :   The now hydrostatic pressure gradient at a given level 
    55       !!      jk is computed by taking the vertical integral of the in-situ  
    56       !!      density gradient along the model level from the suface to that  
    57       !!      level. s-coordinates ('key_s_coord'): a corrective term is added 
    58       !!      to the horizontal pressure gradient : 
    59       !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
    60       !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
     83      !! ** Method  :   Call the hydrostatic pressure gradient routine  
     84      !!      using the scheme defined in the namelist (nhpg parameter) 
     85      !!    
     86      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     87      !!             - Save the trend (l_trddyn=T) 
     88      !!             - Control print  (ln_ctl) 
     89      !! 
     90      !! History : 
     91      !!   9.0  !  05-10  (A. Beckmann, G. Madec) various s-coordinate options 
     92      !!---------------------------------------------------------------------- 
     93      !! * Arguments 
     94      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
     95 
     96      !! * local declarations 
     97      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
     98         ztrdu, ztrdv                         ! 3D temporary workspace 
     99      !!---------------------------------------------------------------------- 
     100    
     101      IF( kt == nit000 )   CALL hpg_ctl      ! initialisation & control of options 
     102 
     103      ! Temporary saving of ua and va trends (l_trddyn) 
     104      IF( l_trddyn )   THEN 
     105         ztrdu(:,:,:) = ua(:,:,:)   
     106         ztrdv(:,:,:) = va(:,:,:)  
     107      ENDIF       
     108 
     109      SELECT CASE ( nhpg )      ! Hydrastatic pressure gradient computation 
     110      CASE ( 0 )                  ! z-coordinate  
     111         CALL hpg_zco( kt ) 
     112      CASE ( 1 )                  ! z-coordinate plus partial steps (interpolation) 
     113         CALL hpg_zps( kt ) 
     114      CASE ( 2 )                  ! s-coordinate (standard jacobian formulation) 
     115         CALL hpg_sco( kt ) 
     116      CASE ( 3 )                  ! s-coordinate (helsinki modification) 
     117         CALL hpg_hel( kt ) 
     118      CASE ( 4 )                  ! s-coordinate (weighted density jacobian) 
     119         CALL hpg_wdj( kt ) 
     120      CASE ( 5 )                  ! s-coordinate (Density Jacobian with Cubic polynomial) 
     121         CALL hpg_djc( kt ) 
     122      CASE ( 6 )                  ! s-coordinate (ROTated axes scheme) 
     123         CALL hpg_rot( kt ) 
     124      CASE ( 10 )                  ! z-coordinate 
     125         CALL hpg_zco_jki( kt ) 
     126      CASE ( 11 )                  ! z-coordinate plus partial steps (interpolation) 
     127         CALL hpg_zps_jki( kt ) 
     128      CASE ( 12 )                  ! s-coordinate (standard jacobian formulation) 
     129         CALL hpg_sco_jki( kt ) 
     130      END SELECT 
     131 
     132      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 
     133      IF( l_trddyn )   THEN 
     134         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     135         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     136         CALL trd_mod( ztrdu, ztrdv, jpdtdhpg, 'DYN', kt ) 
     137      ENDIF           
     138       
     139      IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
     140         CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, & 
     141            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     142      ENDIF 
     143       
     144   END SUBROUTINE dyn_hpg 
     145 
     146 
     147   SUBROUTINE hpg_ctl 
     148      !!---------------------------------------------------------------------- 
     149      !!                 ***  ROUTINE hpg_ctl  *** 
     150      !! 
     151      !! ** Purpose :   initializations for the hydrostatic pressure gradient 
     152      !!              computation and consistency control 
     153      !! 
     154      !! ** Action  :   Read the namelist namdynhpg and check the consistency 
     155      !!      with the type of vertical coordinate used (zco, zps, sco) 
     156      !! 
     157      !! History : 
     158      !!   9.0  !  05-10  (A. Beckmann)  Original code 
     159      !!---------------------------------------------------------------------- 
     160      INTEGER ::   ioptio = 0      ! temporary integer 
     161 
     162      NAMELIST/nam_dynhpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,               & 
     163         &                 ln_hpg_hel, ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot,   & 
     164         &                 gamm 
     165      !!---------------------------------------------------------------------- 
     166 
     167      ! Read Namelist nam_dynhpg : pressure gradient calculation options 
     168      REWIND ( numnam ) 
     169      READ   ( numnam, nam_dynhpg ) 
     170 
     171      ! Control print 
     172      IF(lwp) THEN 
     173         WRITE(numout,*) 
     174         WRITE(numout,*) 'dyn:hpg_ctl : hydrostatic pressure gradient control' 
     175         WRITE(numout,*) '~~~~~~~~~~~' 
     176         WRITE(numout,*) '       Namelist nam_dynhpg : choice of hpg scheme' 
     177         WRITE(numout,*) '          z-coord. - full steps                          ln_hpg_zco = ', ln_hpg_zco 
     178         WRITE(numout,*) '          z-coord. - partial steps (interpolation)       ln_hpg_zps = ', ln_hpg_zps 
     179         WRITE(numout,*) '          s-coord. (standard jacobian formulation)       ln_hpg_sco = ', ln_hpg_sco 
     180         WRITE(numout,*) '          s-coord. (helsinki modification)               ln_hpg_hel = ', ln_hpg_hel 
     181         WRITE(numout,*) '          s-coord. (weighted density jacobian)           ln_hpg_wdj = ', ln_hpg_wdj 
     182         WRITE(numout,*) '          s-coord. (Density Jacobian: Cubic polynomial)  ln_hpg_djc = ', ln_hpg_djc 
     183         WRITE(numout,*) '          s-coord. (ROTated axes scheme)                 ln_hpg_rot = ', ln_hpg_rot 
     184         WRITE(numout,*) '          weighting coeff. (wdj scheme)                     gamm       = ', gamm 
     185      ENDIF 
     186 
     187      ! set nhpg from ln_hpg_... flags 
     188      IF( ln_hpg_zco )   nhpg = 0 
     189      IF( ln_hpg_zps )   nhpg = 1 
     190      IF( ln_hpg_sco )   nhpg = 2 
     191      IF( ln_hpg_hel )   nhpg = 3 
     192      IF( ln_hpg_wdj )   nhpg = 4 
     193      IF( ln_hpg_djc )   nhpg = 5 
     194      IF( ln_hpg_rot )   nhpg = 6 
     195 
     196      ! Consitency check 
     197      ioptio = 0  
     198      IF( ln_hpg_zco )   ioptio = ioptio + 1 
     199      IF( ln_hpg_zps )   ioptio = ioptio + 1 
     200      IF( ln_hpg_sco )   ioptio = ioptio + 1 
     201      IF( ln_hpg_hel )   ioptio = ioptio + 1 
     202      IF( ln_hpg_wdj )   ioptio = ioptio + 1 
     203      IF( ln_hpg_djc )   ioptio = ioptio + 1 
     204      IF( ln_hpg_rot )   ioptio = ioptio + 1 
     205      IF ( ioptio > 1 ) THEN 
     206          IF(lwp) WRITE(numout,cform_err) 
     207          IF(lwp) WRITE(numout,*) ' several hydrostatic pressure gradient options used' 
     208          nstop = nstop + 1 
     209      ENDIF 
     210 
     211      IF( lk_dynhpg_jki ) THEN 
     212         nhpg = nhpg + 10 
     213         IF(lwp) WRITE(numout,*) 
     214         IF(lwp) WRITE(numout,*) '          Autotasking or OPENMP: use j-k-i loops (i.e. _jki routines)' 
     215      ENDIF 
     216 
     217   END SUBROUTINE hpg_ctl 
     218 
     219 
     220   SUBROUTINE hpg_zco( kt ) 
     221      !!--------------------------------------------------------------------- 
     222      !!                  ***  ROUTINE hpg_zco  *** 
     223      !! 
     224      !! ** Method  :   z-coordinate case, levels are horizontal surfaces. 
     225      !!      The now hydrostatic pressure gradient at a given level, jk, 
     226      !!      is computed by taking the vertical integral of the in-situ 
     227      !!      density gradient along the model level from the suface to that 
     228      !!      level:    zhpi = grav ..... 
     229      !!                zhpj = grav ..... 
    61230      !!      add it to the general momentum trend (ua,va). 
    62       !!         ua = ua - 1/e1u * zhpi 
    63       !!         va = va - 1/e2v * zhpj 
    64       !! 
     231      !!            ua = ua - 1/e1u * zhpi 
     232      !!            va = va - 1/e2v * zhpj 
     233      !!  
    65234      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    66       !!             - Save the trend in (utrd,vtrd) ('key_trddyn') 
    67235      !! 
    68236      !! History : 
    69       !!   1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code 
    70       !!        !  91-11  (G. Madec) 
    71       !!        !  96-01  (G. Madec)  s-coordinates 
    72       !!        !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg 
    73       !!   8.5  !  02-08  (G. Madec)  F90: Free form and module, vector opt. 
    74       !!   9.0  !  04-08  (C. Talandier) New trends organization 
     237      !!   1.0  !  87-09  (P. Andrich, M.-A. Foujols)  Original code 
     238      !!   5.0  !  91-11  (G. Madec) 
     239      !!   7.0  !  96-01  (G. Madec) 
     240      !!   8.0  !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg 
     241      !!   8.5  !  02-07  (G. Madec)  F90: Free form and module 
    75242      !!---------------------------------------------------------------------- 
    76243      !! * modules used 
    77244      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace 
    78245         &              zhpj => sa      ! use sa as 3D workspace 
    79  
     246       
    80247      !! * Arguments 
    81248      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
    82249       
    83       !! * Local declarations 
     250      !! * local declarations 
    84251      INTEGER ::   ji, jj, jk           ! dummy loop indices 
    85       REAL(wp) ::   & 
    86          zcoef0, zcoef1, zuap, zvap     ! temporary scalars 
    87       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    88          ztdua, ztdva                   ! temporary scalars 
    89       !!---------------------------------------------------------------------- 
    90  
     252      REAL(wp) ::   &    
     253         zcoef0, zcoef1                 ! temporary scalars 
     254      !!---------------------------------------------------------------------- 
     255       
    91256      IF( kt == nit000 ) THEN 
    92257         IF(lwp) WRITE(numout,*) 
    93          IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend' 
    94          IF(lwp) WRITE(numout,*) '~~~~~~~   s-coordinate case, vector opt. case' 
     258         IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' 
     259         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case ' 
    95260      ENDIF 
    96  
    97       ! Save ua and va trends 
    98       IF( l_trddyn )   THEN 
    99          ztdua(:,:,:) = ua(:,:,:)  
    100          ztdva(:,:,:) = va(:,:,:)  
    101       ENDIF 
    102  
    103       ! 0. Local constant initialization 
    104       ! -------------------------------- 
     261       
     262       
     263      ! Local constant initialization  
     264      ! ----------------------------- 
    105265      zcoef0 = - grav * 0.5 
    106       zuap   = 0.e0 
    107       zvap   = 0.e0 
    108  
    109       ! 1. Surface value 
    110       ! ---------------- 
    111       DO jj = 2, jpjm1 
    112          DO ji = fs_2, fs_jpim1   ! vector opt. 
    113             ! hydrostatic pressure gradient along s-surfaces 
    114             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj)   & 
    115                        * ( fse3w(ji+1,jj,1) * rhd(ji+1,jj,1) - fse3w(ji,jj,1) * rhd(ji,jj,1)  ) 
    116             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj)   & 
    117                        * ( fse3w(ji,jj+1,1) * rhd(ji,jj+1,1) - fse3w(ji,jj,1) * rhd(ji,jj,1)  ) 
    118             ! s-coordinate pressure gradient correction 
    119             zuap = -zcoef0 * ( rhd(ji+1,jj,1) + rhd(ji,jj,1) )   & 
    120                  * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
    121             zvap = -zcoef0 * ( rhd(ji,jj+1,1) + rhd(ji,jj,1) )   & 
    122                  * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
    123             ! add to the general momentum trend 
    124             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    125             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
    126          END DO   
    127       END DO   
    128  
    129       ! 2. interior value (2=<jk=<jpkm1) 
    130       ! ----------------- 
    131       DO jk = 2, jpkm1 
    132          DO jj = 2, jpjm1  
    133             DO ji = fs_2, fs_jpim1   ! vector opt. 
    134                ! hydrostatic pressure gradient along s-surfaces 
    135                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
    136                   &           * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   & 
    137                   &              -fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) 
    138                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
    139                   &           * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   & 
    140                   &              -fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) 
    141                ! s-coordinate pressure gradient correction  
    142                zuap = -zcoef0 * ( rhd(ji+1,jj  ,jk) + rhd(ji,jj,jk) )   & 
    143                     * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
    144                zvap = -zcoef0 * ( rhd(ji  ,jj+1,jk) + rhd(ji,jj,jk) )   & 
    145                     * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
    146                ! add to the general momentum trend 
    147                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    148                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
    149             END DO 
    150          END DO 
    151       END DO 
    152  
    153       ! save the hydrostatic pressure gradient trends for diagnostic 
    154       ! momentum trends 
    155       IF( l_trddyn )   THEN 
    156          zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 
    157          zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:) 
    158          CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt) 
    159       ENDIF 
    160  
    161       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    162          CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, & 
    163             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    164       ENDIF 
    165  
    166    END SUBROUTINE dyn_hpg 
    167  
    168 #elif defined key_partial_steps 
    169    !!--------------------------------------------------------------------- 
    170    !!   'key_partial_steps'                     z-coordinate partial steps 
    171    !!--------------------------------------------------------------------- 
    172  
    173    SUBROUTINE dyn_hpg( kt ) 
    174       !!--------------------------------------------------------------------- 
    175       !!                 ***  ROUTINE dyn_hpg  *** 
    176       !!                     
    177       !! ** Purpose :   Compute the now momentum trend due to the horizontal  
    178       !!      gradient of the hydrostatic pressure. Add it to the general 
    179       !!      momentum trend. 
    180       !! 
    181       !! ** Method  :   The now hydrostatic pressure gradient at a given level  
    182       !!      jk is computed by taking the vertical integral of the in-situ  
    183       !!      density gradient along the model level from the suface to that 
    184       !!      level:   zhpi = grav ..... 
    185       !!               zhpj = grav ..... 
    186       !!      add it to the general momentum trend (ua,va). 
    187       !!            ua = ua - 1/e1u * zhpi 
    188       !!            va = va - 1/e2v * zhpj 
    189       !! 
    190       !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend 
    191       !!              - Save the trend in (utrd,vtrd) ('key_trddyn') 
    192       !! 
    193       !! History : 
    194       !!   8.5  !  02-08  (A. Bozec)  Original code 
    195       !!---------------------------------------------------------------------- 
    196       !! * modules used 
    197       USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace 
    198          &              zhpj => sa      ! use sa as 3D workspace 
    199  
    200       !! * Arguments 
    201       INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
    202  
    203       !! * local declarations 
    204       INTEGER ::   ji, jj, jk           ! dummy loop indices 
    205       INTEGER ::   iku, ikv             ! temporary integers 
    206       REAL(wp) ::   & 
    207          zcoef0, zcoef1, zuap,       &  ! temporary scalars 
    208          zcoef2, zcoef3, zvap           !    "         " 
    209       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    210          ztdua, ztdva                   ! temporary scalars 
    211       !!---------------------------------------------------------------------- 
    212  
    213       IF( kt == nit000 ) THEN 
    214          IF(lwp) WRITE(numout,*) 
    215          IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend' 
    216          IF(lwp) WRITE(numout,*) '~~~~~~~   z-coordinate with partial steps' 
    217          IF(lwp) WRITE(numout,*) '          vector optimization, no autotasking' 
    218       ENDIF 
    219  
    220       ! Save ua and va trends 
    221       IF( l_trddyn )   THEN 
    222          ztdua(:,:,:) = ua(:,:,:)  
    223          ztdva(:,:,:) = va(:,:,:)  
    224       ENDIF 
    225  
    226       ! 0. Local constant initialization 
    227       ! -------------------------------- 
    228       zcoef0 = - grav * 0.5 
    229       zuap   = 0.e0 
    230       zvap   = 0.e0 
    231  
    232       ! 1. Surface value 
    233       ! ---------------- 
     266 
     267      ! Surface value 
     268      ! ------------- 
    234269      DO jj = 2, jpjm1 
    235270         DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    244279      END DO 
    245280 
     281      ! interior value (2=<jk=<jpkm1) 
     282      ! -------------- 
     283      DO jk = 2, jpkm1 
     284         DO jj = 2, jpjm1 
     285            DO ji = fs_2, fs_jpim1   ! vector opt. 
     286               zcoef1 = zcoef0 * fse3w(ji,jj,jk) 
     287               ! hydrostatic pressure gradient 
     288               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
     289                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   & 
     290                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj) 
     291 
     292               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   & 
     293                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   & 
     294                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj) 
     295               ! add to the general momentum trend 
     296               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     297               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
     298            END DO 
     299         END DO 
     300      END DO 
     301 
     302   END SUBROUTINE hpg_zco 
     303 
     304 
     305   SUBROUTINE hpg_zps( kt ) 
     306      !!--------------------------------------------------------------------- 
     307      !!                 ***  ROUTINE hpg_zps  *** 
     308      !!                     
     309      !! ** Method  :   z-coordinate plus partial steps case.  blahblah... 
     310      !!  
     311      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend 
     312      !!     
     313      !! History : 
     314      !!   8.5  !  02-08  (A. Bozec)  Original code 
     315      !!   9.0  !  04-08  (G. Madec)  F90 
     316      !!----------------------------------------------------------------------  
     317      !! * modules used                                  
     318      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace 
     319         &              zhpj => sa      ! use sa as 3D workspace 
     320                
     321      !! * Arguments 
     322      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
     323 
     324      !! * local declarations 
     325      INTEGER ::   ji, jj, jk           ! dummy loop indices 
     326      INTEGER ::   iku, ikv             ! temporary integers 
     327      REAL(wp) ::   & 
     328         zcoef0, zcoef1,             &  ! temporary scalars 
     329         zcoef2, zcoef3                 !    "         " 
     330      !!---------------------------------------------------------------------- 
     331 
     332      IF( kt == nit000 ) THEN 
     333         IF(lwp) WRITE(numout,*) 
     334         IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend' 
     335         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps' 
     336         IF(lwp) WRITE(numout,*) '              vector optimization' 
     337      ENDIF 
     338 
     339 
     340      ! 0. Local constant initialization 
     341      ! -------------------------------- 
     342      zcoef0 = - grav * 0.5 
     343 
     344      ! 1. Surface value 
     345      ! ---------------- 
     346      DO jj = 2, jpjm1 
     347         DO ji = fs_2, fs_jpim1   ! vector opt. 
     348            zcoef1 = zcoef0 * fse3w(ji,jj,1) 
     349            ! hydrostatic pressure gradient 
     350            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 
     351            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) 
     352            ! add to the general momentum trend 
     353            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     354            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     355         END DO 
     356      END DO 
     357 
    246358      ! 2. interior value (2=<jk=<jpkm1) 
    247359      ! ----------------- 
    … …  
    252364               ! hydrostatic pressure gradient 
    253365               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
    254                   &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   & 
    255                   &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj) 
     366                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   & 
     367                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj) 
    256368 
    257369               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   & 
    258                   &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   & 
    259                   &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj) 
     370                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   & 
     371                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj) 
    260372               ! add to the general momentum trend 
    261373               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    262374               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
    263             END DO  
     375            END DO 
    264376         END DO 
    265377      END DO 
    … …  
    279391            ! on i-direction 
    280392            IF ( iku > 2 ) THEN 
    281                ! subtract old value   
     393               ! subtract old value 
    282394               ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) 
    283                ! compute the new one    
     395               ! compute the new one 
    284396               zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1)   & 
    285397                  + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 
    … …  
    289401            ! on j-direction 
    290402            IF ( ikv > 2 ) THEN 
    291                ! subtract old value   
     403               ! subtract old value 
    292404               va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) 
    293                ! compute the new one    
     405               ! compute the new one 
    294406               zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1)   & 
    295407                  + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 
    … …  
    302414      END DO 
    303415 
    304       ! save the hydrostatic pressure gradient trends for diagnostic 
    305       ! momentum trends 
    306       IF( l_trddyn )   THEN 
    307          zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 
    308          zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:) 
    309          CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt) 
    310       ENDIF 
    311  
    312       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    313          CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, & 
    314             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    315       ENDIF 
    316  
    317    END SUBROUTINE dyn_hpg 
    318  
    319 #else 
    320    !!--------------------------------------------------------------------- 
    321    !!   Default case :                                        z-coordinate 
    322    !!--------------------------------------------------------------------- 
    323  
    324    SUBROUTINE dyn_hpg( kt ) 
     416   END SUBROUTINE hpg_zps 
     417 
     418 
     419   SUBROUTINE hpg_sco( kt ) 
    325420      !!--------------------------------------------------------------------- 
    326       !!                  ***  ROUTINE dyn_hpg  *** 
    327       !! 
    328       !! ** Purpose :   Compute the now momentum trend due to the horizontal 
    329       !!      gradient of the hydrostatic pressure. Add it to the general  
    330       !!      momentum trend. 
    331       !! 
    332       !! ** Method  :   The now hydrostatic pressure gradient at a given level 
    333       !!      jk is computed by taking the vertical integral of the in-situ 
     421      !!                  ***  ROUTINE hpg_sco  *** 
     422      !! 
     423      !! ** Method  :   s-coordinate case. Jacobian scheme. 
     424      !!      The now hydrostatic pressure gradient at a given level, jk, 
     425      !!      is computed by taking the vertical integral of the in-situ 
    334426      !!      density gradient along the model level from the suface to that 
    335       !!      level:    zhpi = grav ..... 
    336       !!                zhpj = grav ..... 
     427      !!      level. s-coordinates (ln_sco): a corrective term is added 
     428      !!      to the horizontal pressure gradient : 
     429      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
     430      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
    337431      !!      add it to the general momentum trend (ua,va). 
    338       !!            ua = ua - 1/e1u * zhpi 
    339       !!            va = va - 1/e2v * zhpj 
     432      !!         ua = ua - 1/e1u * zhpi 
     433      !!         va = va - 1/e2v * zhpj 
    340434      !! 
    341435      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    342       !!             - Save the trend in (utrd,vtrd) ('key_trddyn') 
    343436      !! 
    344437      !! History : 
    345       !!   1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code 
    346       !!        !  91-11  (G. Madec) 
    347       !!        !  96-01  (G. Madec)  s-coordinates 
     438      !!   7.0  !  96-01  (G. Madec)  s-coordinates 
    348439      !!        !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg 
    349       !!   8.5  !  02-07  (G. Madec)  F90: Free form and module 
     440      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module, vector opt. 
     441      !!   9.0  !  04-08  (C. Talandier) New trends organization 
     442      !!   9.0  !  05-10  (A. Beckmann) various s-coordinate options 
    350443      !!---------------------------------------------------------------------- 
    351444      !! * modules used 
    … …  
    356449      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
    357450 
    358       !! * local declarations 
     451      !! * Local declarations 
    359452      INTEGER ::   ji, jj, jk           ! dummy loop indices 
    360453      REAL(wp) ::   & 
    361          zcoef0, zcoef1, zuap, zvap     ! temporary scalars 
    362       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    363          ztdua, ztdva                   ! temporary scalars 
     454         zcoef0, zuap, zvap             ! temporary scalars 
    364455      !!---------------------------------------------------------------------- 
    365456 
    366457      IF( kt == nit000 ) THEN 
    367458         IF(lwp) WRITE(numout,*) 
    368          IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend' 
    369          IF(lwp) WRITE(numout,*) '~~~~~~~   z-coordinate case ' 
     459         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 
     460         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
    370461      ENDIF 
    371462 
    372       ! Save ua and va trends 
    373       IF( l_trddyn )   THEN 
    374          ztdua(:,:,:) = ua(:,:,:)  
    375          ztdva(:,:,:) = va(:,:,:)  
    376       ENDIF 
    377463 
    378464      ! 0. Local constant initialization 
    379465      ! -------------------------------- 
    380466      zcoef0 = - grav * 0.5 
    381       zuap   = 0.e0 
    382       zvap   = 0.e0 
    383  
     467 
     468 
     469      ! 1. Surface value 
     470      ! ----------------                           
     471      DO jj = 2, jpjm1 
     472         DO ji = fs_2, fs_jpim1   ! vector opt.    
     473            ! hydrostatic pressure gradient along s-surfaces 
     474            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)   & 
     475               &                                  - fse3w(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) ) 
     476            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   & 
     477               &                                  - fse3w(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) ) 
     478            ! s-coordinate pressure gradient correction 
     479            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   & 
     480               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
     481            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   & 
     482               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     483            ! add to the general momentum trend 
     484            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     485            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     486         END DO   
     487      END DO    
     488             
     489                
     490      ! 2. interior value (2=<jk=<jpkm1) 
     491      ! -----------------      
     492      DO jk = 2, jpkm1                                   
     493         DO jj = 2, jpjm1      
     494            DO ji = fs_2, fs_jpim1   ! vector opt.       
     495               ! hydrostatic pressure gradient along s-surfaces 
     496               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &  
     497                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &  
     498                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) 
     499               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
     500                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   & 
     501                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) 
     502               ! s-coordinate pressure gradient correction 
     503               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) )   & 
     504                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
     505               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) )   & 
     506                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     507               ! add to the general momentum trend 
     508               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     509               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     510            END DO 
     511         END DO 
     512      END DO 
     513 
     514   END SUBROUTINE hpg_sco 
     515 
     516 
     517   SUBROUTINE hpg_hel( kt ) 
     518      !!--------------------------------------------------------------------- 
     519      !!                  ***  ROUTINE hpg_hel  *** 
     520      !! 
     521      !! ** Method  :   s-coordinate case. 
     522      !!      The now hydrostatic pressure gradient at a given level 
     523      !!      jk is computed by taking the vertical integral of the in-situ  
     524      !!      density gradient along the model level from the suface to that  
     525      !!      level. s-coordinates (ln_sco): a corrective term is added 
     526      !!      to the horizontal pressure gradient : 
     527      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
     528      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
     529      !!      add it to the general momentum trend (ua,va). 
     530      !!         ua = ua - 1/e1u * zhpi 
     531      !!         va = va - 1/e2v * zhpj 
     532      !! 
     533      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     534      !!             - Save the trend (l_trddyn=T) 
     535      !! 
     536      !! History : 
     537      !!   9.0  !  05-10  (A. Beckmann)  Original code 
     538      !!---------------------------------------------------------------------- 
     539      !! * modules used 
     540      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace 
     541         &              zhpj => sa      ! use sa as 3D workspace 
     542 
     543      !! * Arguments 
     544      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
     545       
     546      !! * Local declarations 
     547      INTEGER ::   ji, jj, jk           ! dummy loop indices 
     548      REAL(wp) ::   & 
     549         zcoef0, zuap, zvap             ! temporary scalars 
     550      !!---------------------------------------------------------------------- 
     551 
     552      IF( kt == nit000 ) THEN 
     553         IF(lwp) WRITE(numout,*) 
     554         IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend' 
     555         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, helsinki modified scheme' 
     556      ENDIF 
     557 
     558 
     559      ! 0. Local constant initialization 
     560      ! -------------------------------- 
     561      zcoef0 = - grav * 0.5 
     562 
     563  
    384564      ! 1. Surface value 
    385565      ! ---------------- 
    386566      DO jj = 2, jpjm1 
    387567         DO ji = fs_2, fs_jpim1   ! vector opt. 
    388             zcoef1 = zcoef0 * fse3w(ji,jj,1) 
    389             ! hydrostatic pressure gradient 
    390             zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 
    391             zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) 
     568            ! hydrostatic pressure gradient along s-surfaces 
     569            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)  & 
     570               &                                  - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) ) 
     571            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)  & 
     572               &                                  - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) ) 
     573            ! s-coordinate pressure gradient correction 
     574            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   & 
     575               &           * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 
     576            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   & 
     577               &           * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 
    392578            ! add to the general momentum trend 
    393             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
    394             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     579            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     580            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
    395581         END DO 
    396582      END DO 
    … …  
    401587         DO jj = 2, jpjm1 
    402588            DO ji = fs_2, fs_jpim1   ! vector opt. 
    403                zcoef1 = zcoef0 * fse3w(ji,jj,jk) 
    404                ! hydrostatic pressure gradient 
    405                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
    406                   &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   & 
    407                   &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj) 
    408  
    409                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   & 
    410                   &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   & 
    411                   &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj) 
     589               ! hydrostatic pressure gradient along s-surfaces 
     590               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 
     591                  &           +  zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk  ) * rhd(ji+1,jj,jk)     & 
     592                  &                                     -fse3t(ji  ,jj,jk  ) * rhd(ji  ,jj,jk)   ) & 
     593                  &           +  zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1)   & 
     594                  &                                     -fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1) ) 
     595               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 
     596                  &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk)   & 
     597                  &                                     -fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk)   ) & 
     598                  &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) & 
     599                  &                                     -fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1) ) 
     600               ! s-coordinate pressure gradient correction 
     601               zuap = - zcoef0 * ( rhd   (ji+1,jj,jk) + rhd   (ji,jj,jk) )   & 
     602                  &            * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
     603               zvap = - zcoef0 * ( rhd   (ji,jj+1,jk) + rhd   (ji,jj,jk) )   & 
     604                  &            * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 
     605               ! add to the general momentum trend 
     606               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     607               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     608            END DO 
     609         END DO 
     610      END DO 
     611 
     612   END SUBROUTINE hpg_hel 
     613 
     614 
     615   SUBROUTINE hpg_wdj( kt ) 
     616      !!--------------------------------------------------------------------- 
     617      !!                  ***  ROUTINE hpg_wdj  *** 
     618      !! 
     619      !! ** Method  :   Weighted Density Jacobian (wdj) scheme (song 1998) 
     620      !!      The weighting coefficients from the namelist parameter gamm 
     621      !!      (alpha=0.5-gamm ; beta=1-alpha=0.5+gamm) 
     622      !! 
     623      !! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998. 
     624      !! 
     625      !! History : 
     626      !!   9.0  !  05-05  (B.W. An)  Original code 
     627      !!        !  05-10  (G. Madec) style & small optimisation 
     628      !!---------------------------------------------------------------------- 
     629      !! * modules used 
     630      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace 
     631         &              zhpj => sa      ! use sa as 3D workspace 
     632 
     633      !! * Arguments 
     634      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
     635 
     636      !! * Local declarations 
     637      INTEGER ::   ji, jj, jk           ! dummy loop indices 
     638      REAL(wp) ::   & 
     639         zcoef0, zuap, zvap,         &  ! temporary scalars 
     640         zalph , zbeta                  !    "         " 
     641      !!---------------------------------------------------------------------- 
     642 
     643      IF( kt == nit000 ) THEN 
     644         IF(lwp) WRITE(numout,*) 
     645         IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend' 
     646         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   Weighted Density Jacobian' 
     647      ENDIF 
     648 
     649 
     650      ! Local constant initialization 
     651      ! ----------------------------- 
     652      zcoef0 = - grav * 0.5 
     653      zalph  = 0.5 - gamm        ! weighting coefficients (alpha=0.5-gamm) 
     654      zbeta  = 0.5 + gamm        !                        (beta =1-alpha=0.5+gamm) 
     655 
     656      ! Surface value (no ponderation) 
     657      ! ------------- 
     658      DO jj = 2, jpjm1 
     659         DO ji = fs_2, fs_jpim1   ! vector opt. 
     660            ! hydrostatic pressure gradient along s-surfaces 
     661            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * (  fse3w(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)   & 
     662               &                                   - fse3w(ji  ,jj  ,1) * rhd(ji  ,jj  ,1)  ) 
     663            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * (  fse3w(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   & 
     664               &                                   - fse3w(ji  ,jj  ,1) * rhd(ji,  jj  ,1)  ) 
     665            ! s-coordinate pressure gradient correction 
     666            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   & 
     667               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
     668            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   & 
     669               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     670            ! add to the general momentum trend 
     671            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     672            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     673         END DO 
     674      END DO 
     675 
     676      ! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta) 
     677      ! -------------- 
     678      DO jk = 2, jpkm1 
     679         DO jj = 2, jpjm1 
     680            DO ji = fs_2, fs_jpim1   ! vector opt. 
     681               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)                            & 
     682                  &           * (   (            fsde3w(ji+1,jj,jk  ) + fsde3w(ji,jj,jk  )        & 
     683                  &                            - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1)    )   & 
     684                  &               * (  zalph * ( rhd   (ji+1,jj,jk-1) - rhd   (ji,jj,jk-1) )      & 
     685                  &                  + zbeta * ( rhd   (ji+1,jj,jk  ) - rhd   (ji,jj,jk  ) )  )   & 
     686                  &             -   (            rhd   (ji+1,jj,jk  ) + rhd   (ji,jj,jk  )        & 
     687                  &                           - rhd   (ji+1,jj,jk-1) - rhd   (ji,jj,jk-1)     )   & 
     688                  &               * (  zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) )      & 
     689                  &                  + zbeta * ( fsde3w(ji+1,jj,jk  ) - fsde3w(ji,jj,jk  ) )  )  ) 
     690               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)                            & 
     691                  &           * (   (           fsde3w(ji,jj+1,jk  ) + fsde3w(ji,jj,jk  )         & 
     692                  &                           - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1)     )   & 
     693                  &               * (  zalph * ( rhd   (ji,jj+1,jk-1) - rhd   (ji,jj,jk-1) )      & 
     694                  &                  + zbeta * ( rhd   (ji,jj+1,jk  ) - rhd   (ji,jj,jk  ) )  )   & 
     695                  &             -   (            rhd   (ji,jj+1,jk  ) + rhd   (ji,jj,jk  )        & 
     696                  &                            - rhd   (ji,jj+1,jk-1) - rhd   (ji,jj,jk-1)    )   & 
     697                  &               * (  zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) )      & 
     698                  &                  + zbeta * ( fsde3w(ji,jj+1,jk  ) - fsde3w(ji,jj,jk  ) )  )  ) 
    412699               ! add to the general momentum trend 
    413700               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    414701               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
    415             END DO  
    416          END DO 
    417       END DO 
    418  
    419       ! save the hydrostatic pressure ggradient trends for diagnostic 
    420       ! momentum trends 
    421       IF( l_trddyn )   THEN 
    422          zhpi(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 
    423          zhpj(:,:,:) = va(:,:,:) - ztdva(:,:,:) 
    424  
    425          CALL trd_mod(zhpi, zhpj, jpdtdhpg, 'DYN', kt) 
     702            END DO 
     703         END DO 
     704      END DO 
     705 
     706   END SUBROUTINE hpg_wdj 
     707 
     708 
     709   SUBROUTINE hpg_djc( kt ) 
     710      !!--------------------------------------------------------------------- 
     711      !!                  ***  ROUTINE hpg_djc  *** 
     712      !! 
     713      !! ** Method  :   Density Jacobian with Cubic polynomial scheme 
     714      !!  
     715      !! Reference: Shchepetkin, A.F. & J.C. McWilliams, J. Geophys. Res., 
     716      !!            108(C3), 3090, 2003 
     717      !! History : 
     718      !!   9.0  !  05-05  (B.W. An)  Original code 
     719      !!---------------------------------------------------------------------- 
     720      !! * modules used 
     721      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace 
     722         &              zhpj => sa      ! use sa as 3D workspace 
     723 
     724      !! * Arguments 
     725      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
     726       
     727      !! * Local declarations 
     728      INTEGER ::   ji, jj, jk           ! dummy loop indices 
     729      REAL(wp) ::   & 
     730         zcoef0, z1_10, cffu, cffx,  &  ! temporary scalars 
     731                 z1_12, cffv, cffy,  &  !    "         " 
     732         zep   , cffw                   !    "         " 
     733      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &  ! 3D workspace 
     734         drhox, dzx, drhou, dzu, rho_i,     & 
     735         drhoy, dzy, drhov, dzv, rho_j,     & 
     736         drhoz, dzz, drhow, dzw, rho_k 
     737      !!---------------------------------------------------------------------- 
     738 
     739      IF( kt == nit000 ) THEN 
     740         IF(lwp) WRITE(numout,*) 
     741         IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 
     742         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme' 
    426743      ENDIF 
    427744 
    428       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    429          CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, & 
    430             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
     745 
     746      ! 0. Local constant initialization 
     747      ! -------------------------------- 
     748      zcoef0 = - grav * 0.5 
     749      z1_10  = 1.0 / 10.0 
     750      z1_12  = 1.0 / 12.0 
     751 
     752      !---------------------------------------------------------------------------------------- 
     753      !  compute and store in provisional arrays elementary vertical and horizontal differences 
     754      !---------------------------------------------------------------------------------------- 
     755 
     756!!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really 
     757 
     758      DO jk = 2, jpkm1 
     759         DO jj = 2, jpjm1 
     760            DO ji = fs_2, fs_jpim1   ! vector opt. 
     761               drhoz(ji,jj,jk) = rhd   (ji  ,jj  ,jk) - rhd   (ji,jj,jk-1) 
     762               dzz  (ji,jj,jk) = fsde3w(ji  ,jj  ,jk) - fsde3w(ji,jj,jk-1) 
     763               drhox(ji,jj,jk) = rhd   (ji+1,jj  ,jk) - rhd   (ji,jj,jk  ) 
     764               dzx  (ji,jj,jk) = fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk  ) 
     765               drhoy(ji,jj,jk) = rhd   (ji  ,jj+1,jk) - rhd   (ji,jj,jk  ) 
     766               dzy  (ji,jj,jk) = fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk  ) 
     767            END DO 
     768         END DO 
     769      END DO 
     770 
     771      !------------------------------------------------------------------------- 
     772      ! compute harmonic averages using eq. 5.18 
     773      !------------------------------------------------------------------------- 
     774      zep = 1.e-15 
     775 
     776      !!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2) 
     777      !!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj 
     778 
     779      DO jk = 2, jpkm1 
     780         DO jj = 2, jpjm1 
     781            DO ji = fs_2, fs_jpim1   ! vector opt. 
     782               cffw = 2.0 * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1) 
     783 
     784               cffu = 2.0 * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  ) 
     785               cffx = 2.0 * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  ) 
     786   
     787               cffv = 2.0 * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  ) 
     788               cffy = 2.0 * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  ) 
     789 
     790               IF( cffw > zep) THEN 
     791                  drhow(ji,jj,jk) = 2.0 *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   & 
     792                     &                  / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 
     793               ELSE 
     794                  drhow(ji,jj,jk) = 0.e0 
     795               ENDIF 
     796 
     797               dzw(ji,jj,jk) = 2.0 *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   & 
     798                  &                / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 
     799 
     800               IF( cffu > zep ) THEN 
     801                  drhou(ji,jj,jk) = 2.0 *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   & 
     802                     &                  / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 
     803               ELSE 
     804                  drhou(ji,jj,jk ) = 0.e0 
     805               ENDIF 
     806 
     807               IF( cffx > zep ) THEN 
     808                  dzu(ji,jj,jk) = 2.0*dzx(ji+1,jj,jk)*dzx(ji,jj,jk)   & 
     809                     &            /(dzx(ji+1,jj,jk)+dzx(ji,jj,jk)) 
     810               ELSE 
     811                  dzu(ji,jj,jk) = 0.e0 
     812               ENDIF 
     813 
     814               IF( cffv > zep ) THEN 
     815                  drhov(ji,jj,jk) = 2.0 *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   & 
     816                     &                  / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 
     817               ELSE 
     818                  drhov(ji,jj,jk) = 0.e0 
     819               ENDIF 
     820 
     821               IF( cffy > zep ) THEN 
     822                  dzv(ji,jj,jk) = 2.0 *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   & 
     823                     &                / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 
     824               ELSE 
     825                  dzv(ji,jj,jk) = 0.e0 
     826               ENDIF 
     827 
     828            END DO 
     829         END DO 
     830      END DO 
     831 
     832      !---------------------------------------------------------------------------------- 
     833      ! apply boundary conditions at top and bottom using 5.36-5.37 
     834      !---------------------------------------------------------------------------------- 
     835      drhow(:,:, 1 ) = 1.5 * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5 * drhow(:,:,  2  ) 
     836      drhou(:,:, 1 ) = 1.5 * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5 * drhou(:,:,  2  ) 
     837      drhov(:,:, 1 ) = 1.5 * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5 * drhov(:,:,  2  ) 
     838 
     839      drhow(:,:,jpk) = 1.5 * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5 * drhow(:,:,jpkm1) 
     840      drhou(:,:,jpk) = 1.5 * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5 * drhou(:,:,jpkm1) 
     841      drhov(:,:,jpk) = 1.5 * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5 * drhov(:,:,jpkm1) 
     842 
     843 
     844      !-------------------------------------------------------------- 
     845      ! Upper half of top-most grid box, compute and store 
     846      !------------------------------------------------------------- 
     847 
     848!!bug gm   :  e3w-de3w = 0.5*e3w  ....  and de3w(2)-de3w(1)=e3w(2) ....   to be verified 
     849!          true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 
     850 
     851      DO jj = 2, jpjm1 
     852         DO ji = fs_2, fs_jpim1   ! vector opt. 
     853            rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) )            & 
     854               &                   * (  rhd(ji,jj,1)                                 & 
     855               &                     + 0.5 * ( rhd(ji,jj,2) - rhd(ji,jj,1) )         & 
     856               &                           * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )   & 
     857               &                           / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) )  )  
     858         END DO 
     859      END DO 
     860 
     861!!bug gm    : here also, simplification is possible 
     862!!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop 
     863 
     864      DO jk = 2, jpkm1 
     865         DO jj = 2, jpjm1 
     866            DO ji = fs_2, fs_jpim1   ! vector opt. 
     867 
     868               rho_k(ji,jj,jk) = zcoef0 * ( rhd   (ji,jj,jk) + rhd   (ji,jj,jk-1) )                                   & 
     869                  &                     * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) )                                   & 
     870                  &            - grav * z1_10 * (                                                                     & 
     871                  &     ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) )                                                     & 
     872                  &   * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
     873                  &   - ( dzw   (ji,jj,jk) - dzw   (ji,jj,jk-1) )                                                     & 
     874                  &   * ( rhd   (ji,jj,jk) - rhd   (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   & 
     875                  &                             ) 
     876 
     877               rho_i(ji,jj,jk) = zcoef0 * ( rhd   (ji+1,jj,jk) + rhd   (ji,jj,jk) )                                   & 
     878                  &                     * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) )                                   & 
     879                  &            - grav* z1_10 * (                                                                      & 
     880                  &     ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) )                                                     & 
     881                  &   * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
     882                  &   - ( dzu   (ji+1,jj,jk) - dzu   (ji,jj,jk) )                                                     & 
     883                  &   * ( rhd   (ji+1,jj,jk) - rhd   (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   & 
     884                  &                            ) 
     885 
     886               rho_j(ji,jj,jk) = zcoef0 * ( rhd   (ji,jj+1,jk) + rhd   (ji,jj,jk) )                                   & 
     887                  &                     * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) )                                   & 
     888                  &            - grav* z1_10 * (                                                                      & 
     889                  &     ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) )                                                     & 
     890                  &   * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
     891                  &   - ( dzv   (ji,jj+1,jk) - dzv   (ji,jj,jk) )                                                     & 
     892                  &   * ( rhd   (ji,jj+1,jk) - rhd   (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   & 
     893                  &                            ) 
     894 
     895            END DO 
     896         END DO 
     897      END DO 
     898      CALL lbc_lnk(rho_k,'W',1.) 
     899      CALL lbc_lnk(rho_i,'U',1.) 
     900      CALL lbc_lnk(rho_j,'V',1.) 
     901 
     902 
     903      ! --------------- 
     904      !  Surface value 
     905      ! --------------- 
     906      DO jj = 2, jpjm1 
     907         DO ji = fs_2, fs_jpim1   ! vector opt. 
     908            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 
     909            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 
     910            ! add to the general momentum trend 
     911            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     912            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     913         END DO 
     914      END DO 
     915 
     916      ! ---------------- 
     917      !  interior value   (2=<jk=<jpkm1) 
     918      ! ---------------- 
     919      DO jk = 2, jpkm1 
     920         DO jj = 2, jpjm1  
     921            DO ji = fs_2, fs_jpim1   ! vector opt. 
     922               ! hydrostatic pressure gradient along s-surfaces 
     923               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                & 
     924                  &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    & 
     925                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) / e1u(ji,jj) 
     926               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                & 
     927                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
     928                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj) 
     929               ! add to the general momentum trend 
     930               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     931               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
     932            END DO 
     933         END DO 
     934      END DO 
     935 
     936   END SUBROUTINE hpg_djc 
     937 
     938 
     939   SUBROUTINE hpg_rot( kt ) 
     940      !!--------------------------------------------------------------------- 
     941      !!                  ***  ROUTINE hpg_rot  *** 
     942      !! 
     943      !! ** Method  :   rotated axes scheme (Thiem and Berntsen 2005) 
     944      !! 
     945      !! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005. 
     946      !! History : 
     947      !!   9.0  !  05-07 (B.W. An) 
     948      !!   9.0  !  05-10 (A. Beckmann) adapted to non-equidistant and masked grids 
     949      !!---------------------------------------------------------------------- 
     950      !! * modules used 
     951      USE oce, ONLY :   zhpi => ta,         &  ! use ta as 3D workspace 
     952         &              zhpj => sa             ! use sa as 3D workspace 
     953 
     954      !! * Arguments 
     955      INTEGER, INTENT( in ) ::   kt            ! ocean time-step index 
     956       
     957      !! * Local declarations 
     958      INTEGER ::   ji, jj, jk                  ! dummy loop indices 
     959      REAL(wp) ::   & 
     960         zforg, zcoef0, zuap, zmskd1, zmskd1m,                             & 
     961         zfrot        , zvap, zmskd2, zmskd2m 
     962      REAL(wp), DIMENSION(jpi,jpj) ::       &  ! 2D temporary workspace 
     963         zdistr, zsina, zcosa              
     964      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &  ! 3D temporary workspace 
     965         zhpiorg, zhpirot, zhpitra, zhpine, &  
     966         zhpjorg, zhpjrot, zhpjtra, zhpjne 
     967      !!---------------------------------------------------------------------- 
     968 
     969      IF( kt == nit000 ) THEN 
     970         IF(lwp) WRITE(numout,*) 
     971         IF(lwp) WRITE(numout,*) 'dyn:hpg_rot : hydrostatic pressure gradient trend' 
     972         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, rotated axes scheme used' 
    431973      ENDIF 
    432974 
    433    END SUBROUTINE dyn_hpg 
    434  
    435 #endif 
     975      ! ------------------------------- 
     976      !  Local constant initialization 
     977      ! ------------------------------- 
     978      zcoef0 = - grav * 0.5 
     979      zforg  = 0.95e0 
     980      zfrot  = 1.e0 - zforg 
     981 
     982      ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance) 
     983      zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)*e1f(:,:) + e2f(:,:)*e1f(:,:) ) 
     984 
     985      ! sinus and cosinus of diagonal angle at F-point 
     986      zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) ) 
     987      zcosa(:,:) = COS( zsina(:,:) ) 
     988      zsina(:,:) = SIN( zsina(:,:) ) 
     989 
     990      ! --------------- 
     991      !  Surface value 
     992      ! --------------- 
     993      ! compute and add to the general trend the pressure gradients along the axes 
     994      DO jj = 2, jpjm1 
     995         DO ji = fs_2, fs_jpim1   ! vector opt. 
     996            ! hydrostatic pressure gradient along s-surfaces 
     997            zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * (  fse3t(ji+1,jj,1) * rhd(ji+1,jj,1)   & 
     998               &                                      - fse3t(ji  ,jj,1) * rhd(ji  ,jj,1)  ) 
     999            zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * (  fse3t(ji,jj+1,1) * rhd(ji,jj+1,1)   & 
     1000               &                                      - fse3t(ji,jj  ,1) * rhd(ji,jj  ,1)  ) 
     1001            ! s-coordinate pressure gradient correction 
     1002            zuap = -zcoef0 * ( rhd   (ji+1,jj  ,1) + rhd   (ji,jj,1) )   & 
     1003               &           * ( fsdept(ji+1,jj  ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 
     1004            zvap = -zcoef0 * ( rhd   (ji  ,jj+1,1) + rhd   (ji,jj,1) )   & 
     1005               &           * ( fsdept(ji  ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 
     1006            ! add to the general momentum trend 
     1007            ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap ) 
     1008            va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap ) 
     1009         END DO 
     1010      END DO 
     1011 
     1012      ! compute the pressure gradients in the diagonal directions 
     1013      DO jj = 1, jpjm1 
     1014         DO ji = 1, fs_jpim1   ! vector opt. 
     1015            zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji  ,jj,1)      ! mask in the 1st diagnonal 
     1016            zmskd2 = tmask(ji  ,jj+1,1) * tmask(ji+1,jj,1)      ! mask in the 2nd diagnonal 
     1017            ! hydrostatic pressure gradient along s-surfaces 
     1018            zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * (  fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1)   & 
     1019               &                                         - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1)  ) 
     1020            zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * (  fse3t(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   & 
     1021               &                                         - fse3t(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)  ) 
     1022            ! s-coordinate pressure gradient correction 
     1023            zuap = -zdistr(ji,jj) * zmskd1 * ( rhd   (ji+1,jj+1,1) + rhd   (ji  ,jj,1) )   & 
     1024               &                           * ( fsdept(ji+1,jj+1,1) - fsdept(ji  ,jj,1) ) 
     1025            zvap = -zdistr(ji,jj) * zmskd2 * ( rhd   (ji  ,jj+1,1) + rhd   (ji+1,jj,1) )   & 
     1026               &                           * ( fsdept(ji  ,jj+1,1) - fsdept(ji+1,jj,1) ) 
     1027            ! back rotation 
     1028            zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap )   & 
     1029               &            - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 
     1030            zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap )   & 
     1031               &            + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 
     1032         END DO 
     1033      END DO 
     1034 
     1035      ! interpolate and add to the general trend the diagonal gradient 
     1036      DO jj = 2, jpjm1 
     1037         DO ji = fs_2, fs_jpim1   ! vector opt. 
     1038            ! averaging 
     1039            zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji  ,jj-1,1) ) 
     1040            zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj  ,1) ) 
     1041            ! add to the general momentum trend 
     1042            ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1)  
     1043            va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1)  
     1044         END DO 
     1045      END DO 
     1046 
     1047      ! ----------------- 
     1048      ! 2. interior value (2=<jk=<jpkm1) 
     1049      ! ----------------- 
     1050      ! compute and add to the general trend the pressure gradients along the axes 
     1051      DO jk = 2, jpkm1 
     1052         DO jj = 2, jpjm1 
     1053            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1054               ! hydrostatic pressure gradient along s-surfaces 
     1055               zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1)                                                 & 
     1056                  &              +  zcoef0 / e1u(ji,jj) * (  fse3t(ji+1,jj,jk  ) * rhd(ji+1,jj,jk  )   & 
     1057                  &                                        - fse3t(ji  ,jj,jk  ) * rhd(ji  ,jj,jk  )   & 
     1058                  &                                        + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1)   & 
     1059                  &                                        - fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1)  ) 
     1060               zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1)                                                 & 
     1061                  &              +  zcoef0 / e2v(ji,jj) * (  fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk  )   & 
     1062                  &                                        - fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk  )   & 
     1063                  &                                        + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1)   & 
     1064                  &                                        - fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1)  ) 
     1065               ! s-coordinate pressure gradient correction 
     1066               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) )   & 
     1067                  &            * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
     1068               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) )   & 
     1069                  &            * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 
     1070               ! add to the general momentum trend 
     1071               ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap ) 
     1072               va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap ) 
     1073            END DO 
     1074         END DO 
     1075      END DO 
     1076 
     1077      ! compute the pressure gradients in the diagonal directions 
     1078      DO jk = 2, jpkm1 
     1079         DO jj = 1, jpjm1 
     1080            DO ji = 1, fs_jpim1   ! vector opt. 
     1081               zmskd1  = tmask(ji+1,jj+1,jk  ) * tmask(ji  ,jj,jk  )      ! level jk   mask in the 1st diagnonal 
     1082               zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji  ,jj,jk-1)      ! level jk-1    "               "      
     1083               zmskd2  = tmask(ji  ,jj+1,jk  ) * tmask(ji+1,jj,jk  )      ! level jk   mask in the 2nd diagnonal 
     1084               zmskd2m = tmask(ji  ,jj+1,jk-1) * tmask(ji+1,jj,jk-1)      ! level jk-1    "               "      
     1085               ! hydrostatic pressure gradient along s-surfaces 
     1086               zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1)                                                       & 
     1087                  &              + zdistr(ji,jj) * zmskd1  * ( fse3t(ji+1,jj+1,jk  ) * rhd(ji+1,jj+1,jk)     & 
     1088                  &                                           -fse3t(ji  ,jj  ,jk  ) * rhd(ji  ,jj  ,jk) )   & 
     1089                  &              + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1)   & 
     1090                  &                                           -fse3t(ji  ,jj  ,jk-1) * rhd(ji  ,jj  ,jk-1) ) 
     1091               zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1)                                                       & 
     1092                  &              + zdistr(ji,jj) * zmskd2  * ( fse3t(ji  ,jj+1,jk  ) * rhd(ji  ,jj+1,jk)     & 
     1093                  &                                           -fse3t(ji+1,jj  ,jk  ) * rhd(ji+1,jj,  jk) )   & 
     1094                  &              + zdistr(ji,jj) * zmskd2m * ( fse3t(ji  ,jj+1,jk-1) * rhd(ji  ,jj+1,jk-1)   & 
     1095                  &                                           -fse3t(ji+1,jj  ,jk-1) * rhd(ji+1,jj,  jk-1) ) 
     1096               ! s-coordinate pressure gradient correction 
     1097               zuap = - zdistr(ji,jj) * zmskd1 * ( rhd   (ji+1,jj+1,jk) + rhd   (ji  ,jj,jk) )   & 
     1098                  &                            * ( fsdept(ji+1,jj+1,jk) - fsdept(ji  ,jj,jk) ) 
     1099               zvap = - zdistr(ji,jj) * zmskd2 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji+1,jj,jk) )   & 
     1100                  &                            * ( fsdept(ji  ,jj+1,jk) - fsdept(ji+1,jj,jk) ) 
     1101               ! back rotation 
     1102               zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap )   & 
     1103                  &             - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 
     1104               zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap )   & 
     1105                  &             + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 
     1106            END DO 
     1107         END DO 
     1108      END DO 
     1109 
     1110      ! interpolate and add to the general trend 
     1111      DO jk = 2, jpkm1 
     1112         DO jj = 2, jpjm1 
     1113            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1114               ! averaging 
     1115               zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji  ,jj-1,jk) ) 
     1116               zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj  ,jk) ) 
     1117               ! add to the general momentum trend 
     1118               ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk)  
     1119               va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk)  
     1120            END DO 
     1121         END DO 
     1122      END DO 
     1123 
     1124   END SUBROUTINE hpg_rot 
    4361125 
    4371126   !!====================================================================== 

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

    r258 r455  
    1717   USE trdmod_oce      ! ocean variables trends 
    1818   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    19    USE prtctl          ! Print control 
    2019 
    2120   IMPLICIT NONE 
    … …  
    6160      !!              diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ] 
    6261      !!              diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ] 
    63       !!      If lk_sco=F and lk_zps=F, the vertical scale factors in the 
     62      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the 
    6463      !!      rotational part of the diffusion are simplified 
    6564      !!      Add this before trend to the general trend (ua,va): 
    … …  
    7069      !! ** Action : - Update (ua,va) with the before iso-level biharmonic 
    7170      !!               mixing trend. 
    72       !!             - Save in (ztdua,ztdva) the trends ('key_trddyn') 
    7371      !! 
    7472      !! History : 
    … …  
    8179      !!   9.0  !  04-08  (C. Talandier) New trends organization 
    8280      !!---------------------------------------------------------------------- 
    83       !! * Modules used      
    84       USE oce, ONLY :    ztdua => ta,      & ! use ta as 3D workspace    
    85                          ztdva => sa         ! use sa as 3D workspace    
    86  
    8781      !! * Arguments 
    8882      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index 
    … …  
    107101      zlv(:,:) = 0.e0 
    108102 
    109       ! Save ua and va trends 
    110       IF( l_trddyn )   THEN 
    111          ztdua(:,:,:) = ua(:,:,:)  
    112          ztdva(:,:,:) = va(:,:,:)  
    113       ENDIF 
    114103      !                                                ! =============== 
    115104      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    118107         ! --------- 
    119108 
    120          IF( lk_sco .OR. lk_zps ) THEN   ! s-coordinate or z-coordinate with partial steps 
     109         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps 
    121110            zuf(:,:) = rotb(:,:,jk) * fse3f(:,:,jk) 
    122111            DO jj = 2, jpjm1 
    … …  
    129118               END DO 
    130119            END DO 
    131          ELSE                            ! z-coordinate 
     120         ELSE                            ! z-coordinate - full step 
    132121            DO jj = 2, jpjm1 
    133122               DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    162151               zuf(ji,jj) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      & 
    163152                  &                            - zcu(ji  ,jj+1) + zcu(ji,jj)  )   & 
    164 #if defined key_s_coord || defined key_partial_steps 
     153#if defined key_zco 
     154                  &                         / ( e1f(ji,jj)*e2f(ji,jj) ) 
     155#else 
    165156                  &       * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) ) 
    166 #else 
    167                   &                         / ( e1f(ji,jj)*e2f(ji,jj) ) 
    168157#endif 
    169158            END DO   
    … …  
    173162         DO jj = 1, jpjm1 
    174163            DO ji = 1, fs_jpim1   ! vector opt. 
    175 #if defined key_s_coord || defined key_partial_steps 
     164#if defined key_zco 
     165               zlu(ji,jj) = e2u(ji,jj) * zlu(ji,jj) 
     166               zlv(ji,jj) = e1v(ji,jj) * zlv(ji,jj) 
     167#else 
    176168               zlu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj) 
    177169               zlv(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj) 
    178 #else 
    179                zlu(ji,jj) = e2u(ji,jj) * zlu(ji,jj) 
    180                zlv(ji,jj) = e1v(ji,jj) * zlv(ji,jj) 
    181170#endif 
    182171            END DO 
    … …  
    186175         DO jj = 2, jpj 
    187176            DO ji = fs_2, jpi   ! vector opt. 
    188 #if defined key_s_coord || defined key_partial_steps 
     177#if defined key_zco 
     178               zbt = e1t(ji,jj) * e2t(ji,jj) 
     179#else 
    189180               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 
    190 #else 
    191                zbt = e1t(ji,jj) * e2t(ji,jj) 
    192181#endif 
    193182               zut(ji,jj) = (  zlu(ji,jj) - zlu(ji-1,jj  )   & 
    … …  
    207196         DO jj = 2, jpjm1 
    208197            DO ji = fs_2, fs_jpim1   ! vector opt. 
    209 #if defined key_s_coord || defined key_partial_steps 
     198#if defined key_zco 
     199               ze2u = e2u(ji,jj) 
     200               ze2v = e1v(ji,jj) 
     201#else 
    210202               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk) 
    211203               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk) 
    212 #else 
    213                ze2u = e2u(ji,jj) 
    214                ze2v = e1v(ji,jj) 
    215204#endif 
    216205               ! horizontal biharmonic diffusive trends 
    … …  
    229218      END DO                                           !   End of slab 
    230219      !                                                ! =============== 
    231       ! save the lateral diffusion trends for diagnostic 
    232       ! momentum trends 
    233       IF( l_trddyn )   THEN 
    234          ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 
    235          ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) 
    236  
    237          CALL trd_mod(ztdua, ztdva, jpdtdldf, 'DYN', kt) 
    238       ENDIF 
    239  
    240       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    241          CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, & 
    242             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    243       ENDIF 
    244220 
    245221   END SUBROUTINE dyn_ldf_bilap 

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

    r258 r455  
    130130      END DO                                           !   End of slab 
    131131      !                                                ! =============== 
    132  
    133       ! save the lateral diffusion trends for diagnostic 
    134       ! momentum trends 
    135       IF( l_trddyn )   THEN 
    136  
    137          CALL trd_mod(zwk3, zwk4, jpdtdldf, 'DYN', kt) 
    138       ENDIF 
    139  
    140       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    141          CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, & 
    142             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    143       ENDIF 
    144132 
    145133   END SUBROUTINE dyn_ldf_bilapg 

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

    r258 r455  
    2121   USE trdmod_oce      ! ocean variables trends 
    2222   USE ldfslp          ! iso-neutral slopes  
     23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2324   USE in_out_manager  ! I/O manager 
    2425   USE prtctl          ! Print control 
    … …  
    3637   !!---------------------------------------------------------------------- 
    3738   !!   OPA 9.0 , LOCEAN-IPSL (2005)  
    38    !! $Header$  
    39    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    4039   !!---------------------------------------------------------------------- 
    4140 
    … …  
    4645      !!                     ***  ROUTINE dyn_ldf_iso  *** 
    4746      !!                        
    48       !! ** Purpose :   Compute the before trend of the horizontal part of the 
    49       !!      lateral momentum diffusion and add it to the general trend of 
    50       !!      momentum equation. 
     47      !! ** Purpose :   Compute the before trend of the rotated laplacian 
     48      !!      operator of lateral momentum diffusion except the diagonal 
     49      !!      vertical term that will be computed in dynzdf module. Add it 
     50      !!      to the general trend of momentum equation. 
    5151      !! 
    5252      !! ** Method : 
    53       !!        The horizontal component of the lateral diffusive trends on 
    54       !!      momentum is provided by a 2nd order operator rotated along neu- 
    55       !!      tral or geopotential surfaces (in s-coordinates). 
     53      !!        The momentum lateral diffusive trend is provided by a 2nd 
     54      !!      order operator rotated along neutral or geopotential surfaces 
     55      !!      (in s-coordinates). 
    5656      !!      It is computed using before fields (forward in time) and isopyc- 
    57       !!      nal or geopotential slopes computed in routine ldfslp or inildf. 
     57      !!      nal or geopotential slopes computed in routine ldfslp. 
    5858      !!      Here, u and v components are considered as 2 independent scalar 
    5959      !!      fields. Therefore, the property of splitting divergent and rota- 
    … …  
    7676      !!      Add this trend to the general trend (ua,va): 
    7777      !!         ua = ua + diffu 
    78       !!      'key_trddyn' defined: the trends are saved for diagnostics. 
     78      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This 
     79      !!      should be modified for applications others than orca_r2 (!!bug) 
    7980      !! 
    8081      !! ** Action : 
    8182      !!        Update (ua,va) arrays with the before geopotential biharmonic 
    8283      !!      mixing trend. 
    83       !!        Save in (uldftrd,vldftrd) arrays the trends if 'key_trddyn' defined 
     84      !!        Update (avmu,avmv) to accompt for the diagonal vertical component 
     85      !!      of the rotated operator in dynzdf module 
    8486      !! 
    8587      !! History : 
    … …  
    8789      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module 
    8890      !!   9.0  !  04-08  (C. Talandier) New trends organization 
     91      !!        !  05-11  (G. Madec)  s-coordinate: horizontal diffusion 
    8992      !!---------------------------------------------------------------------- 
    90       !! * Modules used      
    91       USE oce, ONLY :    ztdua => ta,   & ! use ta as 3D workspace    
    92                          ztdva => sa      ! use sa as 3D workspace    
    9393      !! * Arguments 
    9494      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index 
    … …  
    103103         ziut, zjuf, zjvt, zivf,        & ! temporary workspace 
    104104         zdku, zdk1u, zdkv, zdk1v 
     105 
     106! ajout dynzdf_iso 
     107      REAL(wp) ::   & 
     108         zcoef0, zcoef3, zcoef4, zmkt, zmkf,   & 
     109         zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj 
     110      REAL(wp), DIMENSION(jpi,jpj) ::        & 
     111         zfuw, zdiu, zdju, zdj1u,            & !    "        " 
     112         zfvw, zdiv, zdjv, zdj1v 
     113 
    105114      !!---------------------------------------------------------------------- 
    106115 
    … …  
    111120      ENDIF 
    112121 
    113       ! Save ua and va trends 
    114       IF( l_trddyn )   THEN 
    115          ztdua(:,:,:) = ua(:,:,:)  
    116          ztdva(:,:,:) = va(:,:,:)  
     122!     ! s-coordinate: Iso-level diffusion on momentum but not on tracer 
     123      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 
     124  
     125         ! set the slopes of iso-level 
     126         DO jk = 1, jpk 
     127            DO jj = 2, jpjm1 
     128               DO ji = fs_2, fs_jpim1   ! vector opt. 
     129                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) 
     130                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
     131                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 
     132                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 
     133               END DO 
     134            END DO 
     135         END DO 
     136  
     137         ! Lateral boundary conditions on the slopes 
     138         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
     139         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
     140  
     141!!bug 
     142         if( kt == nit000 ) then 
     143            IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  & 
     144               &                             ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj)) 
     145         endif 
     146!!end 
    117147      ENDIF 
     148 
    118149      !                                                ! =============== 
    119150      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    142173         ! i-flux at t-point             -----f----- 
    143174 
    144          DO jj = 2, jpjm1 
    145             DO ji = fs_2, jpi   ! vector opt. 
    146                zabe1 = ( fsahmt(ji,jj,jk) + ahmb0 )   & 
    147 #if defined key_partial_steps 
    148                      * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1, jj,jk) ) / e1t(ji,jj) 
    149 #else 
    150                      * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 
    151 #endif 
    152  
    153                zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
    154                               + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
    155  
    156                zcof1 = - aht0 * e2t(ji,jj) * zmskt   & 
    157                      * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    158  
    159                ziut(ji,jj) = tmask(ji,jj,jk) *   & 
    160                            (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
    161                             + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
    162                                        +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) 
    163             END DO 
    164          END DO 
     175         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u) 
     176            DO jj = 2, jpjm1 
     177               DO ji = fs_2, jpi   ! vector opt. 
     178                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj) 
     179 
     180                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
     181                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
     182 
     183                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     184    
     185                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
     186                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
     187                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk) 
     188               END DO 
     189            END DO 
     190         ELSE                   ! other coordinate system (zco or sco) : e3t 
     191            DO jj = 2, jpjm1 
     192               DO ji = fs_2, jpi   ! vector opt. 
     193                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 
     194 
     195                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
     196                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
     197 
     198                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     199 
     200                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
     201                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
     202                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk) 
     203               END DO 
     204            END DO 
     205         ENDIF 
    165206 
    166207         ! j-flux at f-point 
    167208         DO jj = 1, jpjm1 
    168209            DO ji = 1, fs_jpim1   ! vector opt. 
    169                zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 )   & 
    170                      * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 
     210               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 
    171211 
    172212               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   & 
    173                               + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. ) 
    174  
    175                zcof2 = - aht0 * e1f(ji,jj) * zmskf   & 
    176                      * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
    177  
    178                zjuf(ji,jj) = fmask(ji,jj,jk) *   & 
    179                            (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   & 
    180                             + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     & 
    181                                        +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) 
     213                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. ) 
     214 
     215               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
     216 
     217               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   & 
     218                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     & 
     219                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk) 
    182220            END DO 
    183221         END DO 
    … …  
    191229         DO jj = 2, jpjm1 
    192230            DO ji = 1, fs_jpim1   ! vector opt. 
    193                zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 )   & 
    194                      * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
     231               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
    195232 
    196233               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
    197                               + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    198  
    199                zcof1 = - aht0 * e2f(ji,jj) * zmskf   & 
    200                      * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
    201  
    202                zivf(ji,jj) = fmask(ji,jj,jk) *   & 
    203                            (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   & 
    204                             + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     & 
    205                                        +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) 
     234                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ) 
     235 
     236               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
     237 
     238               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   & 
     239                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     & 
     240                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk) 
    206241            END DO 
    207242         END DO 
    208243 
    209244         ! j-flux at t-point 
    210          DO jj = 2, jpj 
    211             DO ji = 1, fs_jpim1   ! vector opt. 
    212                zabe2 = ( fsahmt(ji,jj,jk) + ahmb0 )   & 
    213 #if defined key_partial_steps 
    214                      * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji, jj-1, jk) ) / e2t(ji,jj) 
    215 #else 
    216                      * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 
    217 #endif 
    218  
    219                zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
    220                               + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    221  
    222                zcof2 = - aht0 * e1t(ji,jj) * zmskt   & 
    223                      * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    224  
    225                zjvt(ji,jj) = tmask(ji,jj,jk) *   & 
    226                            (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
    227                             + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     & 
    228                                        +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) 
    229             END DO 
    230          END DO 
     245         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u) 
     246            DO jj = 2, jpj 
     247               DO ji = 1, fs_jpim1   ! vector opt. 
     248                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 
     249 
     250                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
     251                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
     252 
     253                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     254 
     255                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     256                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     & 
     257                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk) 
     258               END DO 
     259            END DO 
     260         ELSE                   ! other coordinate system (zco or sco) : e3t 
     261            DO jj = 2, jpj 
     262               DO ji = 1, fs_jpim1   ! vector opt. 
     263                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 
     264 
     265                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
     266                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
     267 
     268                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     269 
     270                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     271                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     & 
     272                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk) 
     273               END DO 
     274            END DO 
     275         ENDIF 
    231276 
    232277 
    … …  
    241286               ! horizontal component of isopycnal momentum diffusive trends 
    242287               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   & 
    243                        zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu 
     288                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu 
    244289               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   & 
    245                        zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv 
     290                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv 
    246291               ! add the trends to the general trends 
    247292               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 
    … …  
    253298      !                                                ! =============== 
    254299 
    255       ! save the lateral diffusion trends for diagnostic 
    256       ! momentum trends will be saved in dynzdf_iso.F90 
    257       IF( l_trddyn )   THEN 
    258          uldftrd(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 
    259          vldftrd(:,:,:) = va(:,:,:) - ztdva(:,:,:) 
    260       ENDIF 
    261  
    262       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    263          CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, & 
    264             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    265       ENDIF 
     300      ! print sum trends (used for debugging) 
     301      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, & 
     302         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     303 
     304 
     305      !                                                ! =============== 
     306      DO jj = 2, jpjm1                                 !  Vertical slab 
     307         !                                             ! =============== 
     308 
     309  
     310         ! I. vertical trends associated with the lateral mixing 
     311         ! ===================================================== 
     312         !  (excluding the vertical flux proportional to dk[t] 
     313 
     314 
     315         ! I.1 horizontal momentum gradient 
     316         ! -------------------------------- 
     317 
     318         DO jk = 1, jpk 
     319            DO ji = 2, jpi 
     320               ! i-gradient of u at jj 
     321               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) ) 
     322               ! j-gradient of u and v at jj 
     323               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) ) 
     324               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) ) 
     325               ! j-gradient of u and v at jj+1 
     326               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) ) 
     327               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) ) 
     328            END DO 
     329         END DO 
     330         DO jk = 1, jpk 
     331            DO ji = 1, jpim1 
     332               ! i-gradient of v at jj 
     333               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) ) 
     334            END DO 
     335         END DO 
     336 
     337 
     338         ! I.2 Vertical fluxes 
     339         ! ------------------- 
     340 
     341         ! Surface and bottom vertical fluxes set to zero 
     342         DO ji = 1, jpi 
     343            zfuw(ji, 1 ) = 0.e0 
     344            zfvw(ji, 1 ) = 0.e0 
     345            zfuw(ji,jpk) = 0.e0 
     346            zfvw(ji,jpk) = 0.e0 
     347         END DO 
     348 
     349         ! interior (2=<jk=<jpk-1) on U field 
     350         DO jk = 2, jpkm1 
     351            DO ji = 2, jpim1 
     352               zcoef0= 0.5 * aht0 * umask(ji,jj,jk) 
     353 
     354               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 
     355               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 
     356 
     357               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   & 
     358                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. ) 
     359               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   & 
     360                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. ) 
     361 
     362               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi 
     363               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj 
     364               ! vertical flux on u field 
     365               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     & 
     366                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   & 
     367                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     & 
     368                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) 
     369               ! update avmu (add isopycnal vertical coefficient to avmu) 
     370               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 
     371               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0 
     372            END DO 
     373         END DO 
     374 
     375         ! interior (2=<jk=<jpk-1) on V field 
     376         DO jk = 2, jpkm1 
     377            DO ji = 2, jpim1 
     378               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk) 
     379 
     380               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 
     381               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 
     382 
     383               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   & 
     384                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. ) 
     385               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   & 
     386                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. ) 
     387 
     388               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi 
     389               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj 
     390               ! vertical flux on v field 
     391               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     & 
     392                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   & 
     393                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     & 
     394                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) ) 
     395               ! update avmv (add isopycnal vertical coefficient to avmv) 
     396               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 
     397               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0 
     398            END DO 
     399         END DO 
     400 
     401 
     402         ! I.3 Divergence of vertical fluxes added to the general tracer trend 
     403         ! ------------------------------------------------------------------- 
     404         DO jk = 1, jpkm1 
     405            DO ji = 2, jpim1 
     406               ! volume elements 
     407               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
     408               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
     409               ! part of the k-component of isopycnal momentum diffusive trends 
     410               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 
     411               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 
     412               ! add the trends to the general trends 
     413               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 
     414               va(ji,jj,jk) = va(ji,jj,jk) + zvav 
     415            END DO 
     416         END DO 
     417         !                                             ! =============== 
     418      END DO                                           !   End of slab 
     419      !                                                ! =============== 
    266420 
    267421   END SUBROUTINE dyn_ldf_iso 

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

    r258 r455  
    1818   USE trdmod_oce      ! ocean variables trends 
    1919   USE ldfslp          ! iso-neutral slopes  
    20    USE prtctl          ! Print control 
    2120 
    2221   IMPLICIT NONE 
    … …  
    5150      !!         difu = 1/e1u di[ahmt hdivb] - 1/(e2u*e3u) dj-1[e3f ahmf rotb] 
    5251      !!         difv = 1/e2v dj[ahmt hdivb] + 1/(e1v*e3v) di-1[e3f ahmf rotb] 
    53       !!      If 'key_s_coord' key is not activated, the vertical scale factor 
    54       !!      is simplified in the rotational part of the diffusion. 
     52      !!      If lk_zco=T, e3f=e3u=e3v, the vertical scale factor are simplified 
     53      !!      in the rotational part of the diffusion. 
    5554      !!      Add this before trend to the general trend (ua,va): 
    5655      !!            (ua,va) = (ua,va) + (diffu,diffv) 
    … …  
    6059      !! ** Action : - Update (ua,va) with the before iso-level harmonic  
    6160      !!               mixing trend. 
    62       !!             - Save in (ztdua,ztdva) arrays the trends ('key_trddyn') 
    6361      !! 
    6462      !! History : 
    … …  
    6967      !!   9.0  !  04-08 (C. Talandier) New trends organization 
    7068      !!---------------------------------------------------------------------- 
    71       !! * Modules used      
    72       USE oce, ONLY :    ztdua => ta,   & ! use ta as 3D workspace    
    73                          ztdva => sa      ! use sa as 3D workspace    
    74  
    7569      !! * Arguments 
    7670      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index 
    … …  
    8882      ENDIF 
    8983 
    90       ! Save ua and va trends 
    91       IF( l_trddyn )   THEN 
    92          ztdua(:,:,:) = ua(:,:,:)  
    93          ztdva(:,:,:) = va(:,:,:)  
    94       ENDIF 
    95  
    9684      !                                                ! =============== 
    9785      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    9987         DO jj = 2, jpjm1 
    10088            DO ji = fs_2, fs_jpim1   ! vector opt. 
    101 #if defined key_s_coord || defined key_partial_steps 
     89#if defined key_zco 
     90               ! horizontal diffusive trends 
     91               ze2u = rotb (ji,jj,jk)*fsahmf(ji,jj,jk) 
     92               ze1v = hdivb(ji,jj,jk)*fsahmt(ji,jj,jk) 
     93               zua = - ( ze2u - rotb (ji,jj-1,jk)*fsahmf(ji,jj-1,jk)                   ) / e2u(ji,jj)   & 
     94                     + ( hdivb(ji+1,jj,jk)*fsahmt(ji+1,jj,jk) - ze1v                   ) / e1u(ji,jj) 
     95 
     96               zva = + ( ze2u - rotb (ji-1,jj,jk)*fsahmf(ji-1,jj,jk)                   ) / e1v(ji,jj)   & 
     97                     + ( hdivb(ji,jj+1,jk)*fsahmt(ji,jj+1,jk) - ze1v                   ) / e2v(ji,jj) 
     98#else 
    10299               ze2u = rotb (ji,jj,jk)*fsahmf(ji,jj,jk)*fse3f(ji,jj,jk) 
    103100               ze1v = hdivb(ji,jj,jk)*fsahmt(ji,jj,jk) 
    … …  
    108105               zva = + ( ze2u - rotb (ji-1,jj,jk)*fsahmf(ji-1,jj,jk)*fse3f(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   & 
    109106                     + ( hdivb(ji,jj+1,jk)*fsahmt(ji,jj+1,jk) - ze1v                   ) / e2v(ji,jj) 
    110 #else 
    111                ! horizontal diffusive trends 
    112                ze2u = rotb (ji,jj,jk)*fsahmf(ji,jj,jk) 
    113                ze1v = hdivb(ji,jj,jk)*fsahmt(ji,jj,jk) 
    114                zua = - (                ze2u                  - rotb (ji,jj-1,jk)*fsahmf(ji,jj-1,jk) ) / e2u(ji,jj)   & 
    115                      + ( hdivb(ji+1,jj,jk)*fsahmt(ji+1,jj,jk) -                ze1v                  ) / e1u(ji,jj) 
    116  
    117                zva = + (                ze2u                  - rotb (ji-1,jj,jk)*fsahmf(ji-1,jj,jk) ) / e1v(ji,jj)   & 
    118                      + ( hdivb(ji,jj+1,jk)*fsahmt(ji,jj+1,jk) -                ze1v                  ) / e2v(ji,jj) 
    119107#endif 
    120108 
    … …  
    128116      !                                                ! =============== 
    129117 
    130       ! save the lateral diffusion trends for diagnostic 
    131       ! momentum trends 
    132       IF( l_trddyn )   THEN 
    133          ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 
    134          ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) 
    135  
    136          CALL trd_mod(ztdua, ztdva, jpdtdldf, 'DYN', kt) 
    137       ENDIF 
    138  
    139       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    140          CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, & 
    141             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    142       ENDIF 
    143  
    144118   END SUBROUTINE dyn_ldf_lap 
    145119 

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

    r371 r455  
    44   !! Ocean dynamics:  surface pressure gradient trend 
    55   !!====================================================================== 
    6 #if ( defined key_dynspg_exp && ! defined key_autotasking ) ||   defined key_esopa 
    7    !!---------------------------------------------------------------------- 
    8    !!   'key_dynspg_exp'       free sfce cst vol. without filter nor ts 
    9    !!   NOT 'key_autotasking'                      k-j-i loop (vector opt.) 
     6#if ( defined key_dynspg_exp && ! defined key_mpp_omp ) ||   defined key_esopa 
     7   !!---------------------------------------------------------------------- 
     8   !!   'key_dynspg_exp'           free sfce cst vol. without filter nor ts 
     9   !!   NOT 'key_mpp_omp'                          k-j-i loop (vector opt.) 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   dyn_spg_exp  : update the momentum trend with the surface  

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

    r371 r455  
    44   !! Ocean dynamics:  surface pressure gradient trend 
    55   !!====================================================================== 
    6 #if ( defined key_dynspg_exp && defined key_autotasking )   ||   defined key_esopa 
     6#if ( defined key_dynspg_exp && defined key_mpp_omp )   ||   defined key_esopa 
    77   !!---------------------------------------------------------------------- 
    88   !!   'key_dynspg_exp'                              explicit free surface 
    9    !!   'key_autotasking'                                        j-k-i loop 
     9   !!   'key_mpp_omp'                                            j-k-i loop 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   dyn_spg_exp_jki  : update the momentum trend with the surface  

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

    r413 r455  
    44   !! Ocean dynamics:  surface pressure gradient trend 
    55   !!====================================================================== 
    6 #if ( defined key_dynspg_flt && ! defined key_autotasking )   ||   defined key_esopa 
     6#if ( defined key_dynspg_flt && ! defined key_mpp_omp )   ||   defined key_esopa 
    77   !!---------------------------------------------------------------------- 
    88   !!   'key_dynspg_flt'                              filtered free surface 
    9    !!   NOT 'key_autotasking'                      k-j-i loop (vector opt.) 
     9   !!   NOT 'key_mpp_omp'                          k-j-i loop (vector opt.) 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   dyn_spg_flt  : update the momentum trend with the surface pressure 

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

    r413 r455  
    44   !! Ocean dynamics:  surface pressure gradient trend 
    55   !!====================================================================== 
    6 #if ( defined key_dynspg_flt && defined key_autotasking )   ||   defined key_esopa 
    7    !!---------------------------------------------------------------------- 
    8    !!   'key_dynspg_flt' & 'key_autotasking'          filtered free surface 
    9    !!                                                   j-k-i loop (j-slab) 
     6#if ( defined key_dynspg_flt && defined key_mpp_omp )   ||   defined key_esopa 
     7   !!---------------------------------------------------------------------- 
     8   !!   'key_dynspg_flt'                              filtered free surface 
     9   !!   'key_mpp_omp'                              j-k-i loop (vector opt.) 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   dyn_spg_flt_jki : Update the momentum trend with the surface pressure 

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

    r374 r455  
    44   !! Ocean dynamics:  surface pressure gradient trend 
    55   !!====================================================================== 
    6 #if ( defined key_dynspg_ts && ! defined key_autotasking ) ||   defined key_esopa 
    7    !!---------------------------------------------------------------------- 
    8    !!   'key_dynspg_ts'     free surface cst volume with time splitting 
    9    !!   NOT 'key_autotasking'                      k-j-i loop (vector opt.) 
     6#if ( defined key_dynspg_ts && ! defined key_mpp_omp ) ||   defined key_esopa 
     7   !!---------------------------------------------------------------------- 
     8   !!   'key_dynspg_ts'         free surface cst volume with time splitting 
     9   !!   NOT 'key_mpp_omp'                          k-j-i loop (vector opt.) 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 

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

    r374 r455  
    44   !! Ocean dynamics:  surface pressure gradient trend 
    55   !!====================================================================== 
    6 #if ( defined key_dynspg_ts && defined key_autotasking )   ||   defined key_esopa 
     6#if ( defined key_dynspg_ts && defined key_mpp_omp )   ||   defined key_esopa 
    77   !!---------------------------------------------------------------------- 
    88   !!   'key_dynspg_ts'                    free surface with time splitting 
    9    !!   'key_autotasking'                          j-k-i loop (vector opt.) 
     9   !!   'key_mpp_omp'                              j-k-i loop (vector opt.) 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 

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

    r418 r455  
    77 
    88   !!---------------------------------------------------------------------- 
    9    !!   dyn_vor_enstrophy: enstrophy conserving scheme       (ln_dynvor_ens=T) 
    10    !!   dyn_vor_energy   : energy conserving scheme          (ln_dynvor_ene=T) 
    11    !!   dyn_vor_mixed    : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 
    12    !!   dyn_vor_ene_ens  : energy and enstrophy conserving   (ln_dynvor_een=T) 
    13    !!   dyn_vor_ctl      : control of the different vorticity option 
     9   !!   dyn_vor     : Update the momentum trend with the vorticity trend 
     10   !!       vor_ens : enstrophy conserving scheme       (ln_dynvor_ens=T) 
     11   !!       vor_ene : energy conserving scheme          (ln_dynvor_ene=T) 
     12   !!       vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 
     13   !!       vor_een : energy and enstrophy conserving   (ln_dynvor_een=T) 
     14   !!       vor_ctl : control of the different vorticity option 
    1415   !!---------------------------------------------------------------------- 
    1516   !! * Modules used 
    … …  
    2627 
    2728   !! * Routine accessibility 
    28    PUBLIC dyn_vor_enstrophy      ! routine called by step.F90 
    29    PUBLIC dyn_vor_energy         ! routine called by step.F90 
    30    PUBLIC dyn_vor_mixed          ! routine called by step.F90 
    31    PUBLIC dyn_vor_ene_ens        ! routine called by step.F90 
    32    PUBLIC dyn_vor_ctl            ! routine called by step.F90 
     29   PUBLIC dyn_vor      ! routine called by step.F90 
    3330 
    3431   !! * Shared module variables 
    … …  
    3835   LOGICAL, PUBLIC ::   ln_dynvor_een = .FALSE.    !: energy and enstrophy conserving scheme 
    3936 
     37   !! * module variables 
     38   INTEGER ::                        & 
     39      nvor = 0                         ! type of vorticity trend used 
     40 
    4041   !! * Substitutions 
    4142#  include "domzgr_substitute.h90" 
    … …  
    4344   !!---------------------------------------------------------------------- 
    4445   !!   OPA 9.0 , LOCEAN-IPSL (2005)  
    45    !! $Header$  
    46    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    4746   !!---------------------------------------------------------------------- 
    4847 
    4948CONTAINS 
    5049 
    51    SUBROUTINE dyn_vor_energy( kt ) 
    52       !!---------------------------------------------------------------------- 
    53       !!                  ***  ROUTINE dyn_vor_energy  *** 
     50   SUBROUTINE dyn_vor( kt ) 
     51      !!---------------------------------------------------------------------- 
     52      !! 
     53      !! ** Purpose :   compute the lateral ocean tracer physics. 
     54      !! 
     55      !! ** Action : - Update (ua,va) with the now vorticity term trend 
     56      !!             - save the trends in (utrd,vtrd) in 2 parts (relative 
     57      !!               and planetary vorticity trends) ('key_trddyn') 
     58      !! 
     59      !! History : 
     60      !!   9.0  !  05-11  (G. Madec)  Original code 
     61      !!---------------------------------------------------------------------- 
     62      USE oce, ONLY :    ztrdu => ta,     & ! use ta as 3D workspace 
     63                         ztrdv => sa        ! use sa as 3D workspace 
     64 
     65      !! * Arguments 
     66      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
     67      !!---------------------------------------------------------------------- 
     68 
     69      IF( kt == nit000 )   CALL vor_ctl          ! initialisation & control of options 
     70 
     71      !                                          ! vorticity term including Coriolis 
     72      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend 
     73 
     74      CASE ( -1 )                                      ! esopa: test all possibility with control print 
     75         CALL vor_ene( kt, 'TOT', ua, va ) 
     76         IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 
     77            &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     78         CALL vor_ens( kt, 'TOT', ua, va ) 
     79         IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 
     80            &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     81         CALL vor_mix( kt ) 
     82         IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 
     83            &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     84         CALL vor_een( kt, 'TOT', ua, va ) 
     85         IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 
     86            &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     87 
     88      CASE ( 0 )                                       ! energy conserving scheme 
     89         IF( l_trddyn )   THEN 
     90            ztrdu(:,:,:) = ua(:,:,:) 
     91            ztrdv(:,:,:) = va(:,:,:) 
     92            CALL vor_ene( kt, 'VOR', ua, va )                ! relative vorticity trend 
     93            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     94            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     95            CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 
     96            ztrdu(:,:,:) = ua(:,:,:) 
     97            ztrdv(:,:,:) = va(:,:,:) 
     98            CALL vor_ene( kt, 'COR', ua, va )                ! planetary vorticity trend 
     99            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     100            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     101            CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 
     102            CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 
     103         ELSE 
     104            CALL vor_ene( kt, 'TOT', ua, va )                ! total vorticity 
     105         ENDIF 
     106 
     107      CASE ( 1 )                                       ! enstrophy conserving scheme 
     108         IF( l_trddyn )   THEN     
     109            ztrdu(:,:,:) = ua(:,:,:) 
     110            ztrdv(:,:,:) = va(:,:,:) 
     111            CALL vor_ens( kt, 'VOR', ua, va )                ! relative vorticity trend 
     112            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     113            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     114            CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 
     115            ztrdu(:,:,:) = ua(:,:,:) 
     116            ztrdv(:,:,:) = va(:,:,:) 
     117            CALL vor_ens( kt, 'COR', ua, va )                ! planetary vorticity trend 
     118            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     119            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     120            CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 
     121            CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 
     122         ELSE 
     123            CALL vor_ens( kt, 'TOT', ua, va )                ! total vorticity 
     124         ENDIF 
     125 
     126      CASE ( 2 )                                       ! mixed ene-ens scheme 
     127         IF( l_trddyn )   THEN 
     128            ztrdu(:,:,:) = ua(:,:,:) 
     129            ztrdv(:,:,:) = va(:,:,:) 
     130            CALL vor_ens( kt, 'VOR', ua, va )                ! relative vorticity trend (ens) 
     131            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     132            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     133            CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 
     134            ztrdu(:,:,:) = ua(:,:,:) 
     135            ztrdv(:,:,:) = va(:,:,:) 
     136            CALL vor_ene( kt, 'COR', ua, va )                ! planetary vorticity trend (ene) 
     137            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     138            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     139            CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 
     140            CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 
     141         ELSE 
     142            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene) 
     143         ENDIF 
     144 
     145      CASE ( 3 )                                       ! energy and enstrophy conserving scheme 
     146         IF( l_trddyn )   THEN 
     147            ztrdu(:,:,:) = ua(:,:,:) 
     148            ztrdv(:,:,:) = va(:,:,:) 
     149            CALL vor_een( kt, 'VOR', ua, va )                ! relative vorticity trend 
     150            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     151            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     152            CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 
     153            ztrdu(:,:,:) = ua(:,:,:) 
     154            ztrdv(:,:,:) = va(:,:,:) 
     155            CALL vor_een( kt, 'COR', ua, va )                ! planetary vorticity trend 
     156            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     157            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     158            CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 
     159            CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 
     160         ELSE 
     161            CALL vor_een( kt, 'TOT', ua, va )                ! total vorticity 
     162         ENDIF 
     163 
     164      END SELECT 
     165 
     166      !                       ! print sum trends (used for debugging) 
     167      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, & 
     168         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     169 
     170   END SUBROUTINE dyn_vor 
     171 
     172 
     173   SUBROUTINE vor_ene( kt, cd_vor, pua, pva ) 
     174      !!---------------------------------------------------------------------- 
     175      !!                  ***  ROUTINE vor_ene  *** 
    54176      !! 
    55177      !! ** Purpose :   Compute the now total vorticity trend and add it to  
    … …  
    60182      !!      horizontal kinetic energy. 
    61183      !!      The trend of the vorticity term is given by: 
    62       !!       * s-coordinate (lk_sco=T), the e3. are inside the derivatives: 
     184      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 
    63185      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ] 
    64186      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ] 
    … …  
    81203      !!   9.0  !  04-08  (C. Talandier)  New trends organization 
    82204      !!---------------------------------------------------------------------- 
    83       !! * Modules used      
    84       USE oce, ONLY :    ztdua => ta,     & ! use ta as 3D workspace    
    85                          ztdva => sa        ! use sa as 3D workspace    
    86  
    87205      !! * Arguments 
    88206      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index 
     207      CHARACTER(len=3) , INTENT( in ) ::   & 
     208         cd_vor        ! define the vorticity considered 
     209         !             ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 
     210      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   & 
     211         pua, pva                      ! ???? 
    89212 
    90213      !! * Local declarations 
    91214      INTEGER ::   ji, jj, jk               ! dummy loop indices 
    92215      REAL(wp) ::   & 
    93          zfact2, zua, zva,               &  ! temporary scalars 
     216         zfact2,                         &  ! temporary scalars 
    94217         zx1, zx2, zy1, zy2                 !    "         " 
    95218      REAL(wp), DIMENSION(jpi,jpj) ::   & 
    96219         zwx, zwy, zwz                      ! temporary workspace 
    97       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    98          zcu, zcv                           !    "         " 
    99220      !!---------------------------------------------------------------------- 
    100221 
    101222      IF( kt == nit000 ) THEN 
    102223         IF(lwp) WRITE(numout,*) 
    103          IF(lwp) WRITE(numout,*) 'dyn_vor_energy : vorticity term: energy conserving scheme' 
    104          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 
     224         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' 
     225         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    105226      ENDIF 
    106227 
    107228      ! Local constant initialization 
    108229      zfact2 = 0.5 * 0.5 
    109       zcu (:,:,:) = 0.e0 
    110       zcv (:,:,:) = 0.e0 
    111  
    112       ! Save ua and va trends 
    113       IF( l_trddyn )   THEN 
    114          ztdua(:,:,:) = ua(:,:,:)  
    115          ztdva(:,:,:) = va(:,:,:)  
    116          zcu(:,:,:) = 0.e0 
    117          zcv(:,:,:) = 0.e0 
    118       ENDIF 
    119        
     230 
     231!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
     232!$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
    120233      !                                                ! =============== 
    121234      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    124237         ! Potential vorticity and horizontal fluxes 
    125238         ! ----------------------------------------- 
    126          IF( lk_sco ) THEN 
    127             zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) / fse3f(:,:,jk) 
     239         SELECT CASE( cd_vor )      ! vorticity considered 
     240         CASE ( 'COR' )                   ! planetary vorticcity (Coriolis) 
     241            zwz(:,:) = ff(:,:) 
     242         CASE ( 'VOR' )                   ! relative  vorticity 
     243            zwz(:,:) = rotn(:,:,jk) 
     244         CASE ( 'TOT' )                   ! total     vorticity (planetary + relative) 
     245            zwz(:,:) = rotn(:,:,jk) + ff(:,:) 
     246         END SELECT 
     247 
     248         IF( ln_sco ) THEN 
     249            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 
    128250            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    129251            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    130252         ELSE 
    131             zwz(:,:) = rotn(:,:,jk) + ff(:,:) 
    132253            zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    133254            zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
    … …  
    142263               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    143264               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    144                zua = zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    145                zva =-zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    146                ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    147                va(ji,jj,jk) = va(ji,jj,jk) + zva 
     265               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     266               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
    148267            END DO   
    149268         END DO   
    … …  
    152271      !                                                ! =============== 
    153272 
    154       ! save the relative & planetary vorticity trends for diagnostic 
    155       ! momentum trends 
    156       IF( l_trddyn )   THEN 
    157          ! Compute the planetary vorticity term trend 
    158          !                                                ! =============== 
    159          DO jk = 1, jpkm1                                 ! Horizontal slab 
    160             !                                             ! =============== 
    161             DO jj = 2, jpjm1 
    162                DO ji = fs_2, fs_jpim1   ! vector opt. 
    163                   zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
    164                   zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
    165                   zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    166                   zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    167 # if defined key_s_coord 
    168                  zcu(ji,jj,jk) = zfact2 / e1u(ji,jj) * ( ff(ji  ,jj-1) / fse3f(ji,jj-1,jk) * zy1  & 
    169                    &                                   + ff(ji  ,jj  ) / fse3f(ji,jj  ,jk) * zy2 ) 
    170                  zcv(ji,jj,jk) =-zfact2 / e2v(ji,jj) * ( ff(ji-1,jj  ) / fse3f(ji-1,jj,jk) * zx1  & 
    171                    &                                   + ff(ji  ,jj  ) / fse3f(ji  ,jj,jk) * zx2 ) 
    172 # else 
    173                  zcu(ji,jj,jk) = zfact2 / e1u(ji,jj) * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 
    174                  zcv(ji,jj,jk) =-zfact2 / e2v(ji,jj) * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) 
    175 # endif 
    176                END DO   
    177             END DO   
    178             !                                             ! =============== 
    179          END DO                                           !   End of slab 
    180          !                                                ! =============== 
    181  
    182          ! Compute the relative vorticity term trend 
    183          ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:)  
    184          ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:)  
    185  
    186          CALL trd_mod(zcu  , zcv  , jpdtdpvo, 'DYN', kt) 
    187          CALL trd_mod(zcu  , zcv  , jpdtddat, 'DYN', kt) 
    188          CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 
    189       ENDIF 
    190  
    191       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    192          CALL prt_ctl(tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, & 
    193             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    194       ENDIF 
    195  
    196    END SUBROUTINE dyn_vor_energy 
    197  
    198  
    199    SUBROUTINE dyn_vor_mixed( kt ) 
    200       !!---------------------------------------------------------------------- 
    201       !!                 ***  ROUTINE dyn_vor_mixed  *** 
     273   END SUBROUTINE vor_ene 
     274 
     275 
     276   SUBROUTINE vor_mix( kt ) 
     277      !!---------------------------------------------------------------------- 
     278      !!                 ***  ROUTINE vor_mix  *** 
    202279      !! 
    203280      !! ** Purpose :   Compute the now total vorticity trend and add it to 
    … …  
    209286      !!      ticity term and the horizontal kinetic energy for (f x uh), the 
    210287      !!      coriolis term. the now trend of the vorticity term is given by: 
    211       !!       * s-coordinate (lk_sco=T), the e3. are inside the derivatives: 
     288      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 
    212289      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ] 
    213290      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ] 
    … …  
    234311      !!   9.0  !  04-08  (C. Talandier)  New trends organization 
    235312      !!---------------------------------------------------------------------- 
    236       !! * Modules used      
    237       USE oce, ONLY :    ztdua => ta,     & ! use ta as 3D workspace    
    238                          ztdva => sa        ! use sa as 3D workspace    
    239313      !! * Arguments 
    240314      INTEGER, INTENT( in ) ::   kt         ! ocean timestep index 
    … …  
    247321      REAL(wp), DIMENSION(jpi,jpj) ::   & 
    248322         zwx, zwy, zwz, zww                 ! temporary workspace 
    249       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    250          zcu, zcv                           !    "         " 
    251323      !!---------------------------------------------------------------------- 
    252324 
    253325      IF( kt == nit000 ) THEN 
    254326         IF(lwp) WRITE(numout,*) 
    255          IF(lwp) WRITE(numout,*) 'dyn_vor_mixed : vorticity term: mixed energy/enstrophy conserving scheme' 
    256          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 
     327         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme' 
     328         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    257329      ENDIF 
    258330 
    … …  
    261333      zfact2 = 0.5 * 0.5 
    262334 
    263       ! Save ua and va trends 
    264       IF( l_trddyn )   THEN 
    265          ztdua(:,:,:) = ua(:,:,:)  
    266          ztdva(:,:,:) = va(:,:,:)  
    267          zcu(:,:,:) = 0.e0 
    268          zcv(:,:,:) = 0.e0 
    269       ENDIF 
    270  
     335!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) 
     336!$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) 
    271337      !                                                ! =============== 
    272338      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    275341         ! Relative and planetary potential vorticity and horizontal fluxes 
    276342         ! ---------------------------------------------------------------- 
    277          IF( lk_sco ) THEN         
     343         IF( ln_sco ) THEN         
    278344            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk) 
    279345            zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk) 
    … …  
    310376      !                                                ! =============== 
    311377 
    312       ! save the relative & planetary vorticity trends for diagnostic 
    313       ! momentum trends 
    314       IF( l_trddyn )   THEN 
    315          ! Compute the planetary vorticity term trend 
    316          !                                                ! =============== 
    317          DO jk = 1, jpkm1                                 ! Horizontal slab 
    318             !                                             ! =============== 
    319             DO jj = 2, jpjm1 
    320                DO ji = fs_2, fs_jpim1   ! vector opt. 
    321                   zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 
    322                   zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    323                   zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 
    324                   zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    325                   ! energy conserving formulation for planetary vorticity term 
    326                   zcu(ji,jj,jk) = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    327                   zcv(ji,jj,jk) =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    328                END DO   
    329             END DO   
    330             !                                             ! =============== 
    331          END DO                                           !   End of slab 
    332          !                                                ! =============== 
    333  
    334          ! Compute the relative vorticity term trend 
    335          ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:)  
    336          ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:)  
    337  
    338          CALL trd_mod(zcu  , zcv  , jpdtdpvo, 'DYN', kt) 
    339          CALL trd_mod(zcu  , zcv  , jpdtddat, 'DYN', kt) 
    340          CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 
    341       ENDIF 
    342  
    343       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    344          CALL prt_ctl(tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, & 
    345             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    346       ENDIF 
    347  
    348    END SUBROUTINE dyn_vor_mixed 
    349  
    350  
    351    SUBROUTINE dyn_vor_enstrophy( kt ) 
    352       !!---------------------------------------------------------------------- 
    353       !!                ***  ROUTINE dyn_vor_enstrophy  *** 
     378   END SUBROUTINE vor_mix 
     379 
     380 
     381   SUBROUTINE vor_ens( kt, cd_vor, pua, pva ) 
     382      !!---------------------------------------------------------------------- 
     383      !!                ***  ROUTINE vor_ens  *** 
    354384      !! 
    355385      !! ** Purpose :   Compute the now total vorticity trend and add it to 
    … …  
    360390      !!      potential enstrophy of a horizontally non-divergent flow. the 
    361391      !!      trend of the vorticity term is given by: 
    362       !!       * s-coordinate (lk_sco=T), the e3. are inside the derivative: 
     392      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative: 
    363393      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ] 
    364394      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ] 
    … …  
    381411      !!   9.0  !  04-08  (C. Talandier)  New trends organization 
    382412      !!---------------------------------------------------------------------- 
    383       !! * modules used 
    384       USE oce, ONLY:   zwx  => ta,   & ! use ta as 3D workspace 
    385                        zwy  => sa      ! use sa as 3D workspace 
    386413      !! * Arguments 
    387414      INTEGER, INTENT( in ) ::   kt    ! ocean timestep 
     415      CHARACTER(len=3) , INTENT( in ) ::   & 
     416         cd_vor        ! define the vorticity considered 
     417         !             ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 
     418      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   & 
     419         pua, pva                      ! ???? 
    388420 
    389421      !! * Local declarations 
    390422      INTEGER ::   ji, jj, jk          ! dummy loop indices 
    391423      REAL(wp) ::   & 
    392          zfact1, zua, zva, zuav, zvau  ! temporary scalars 
    393       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    394          zcu, zcv, zwz,              & ! temporary workspace 
    395          ztdua, ztdva                  ! temporary workspace 
     424         zfact1, zuav, zvau            ! temporary scalars 
     425      REAL(wp), DIMENSION(jpi,jpj) ::   & 
     426         zwx, zwy, zwz                 ! temporary workspace 
    396427      !!---------------------------------------------------------------------- 
    397428       
    398429      IF( kt == nit000 ) THEN 
    399430         IF(lwp) WRITE(numout,*) 
    400          IF(lwp) WRITE(numout,*) 'dyn_vor_enstrophy : vorticity term: enstrophy conserving scheme' 
    401          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 
     431         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme' 
     432         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    402433      ENDIF 
    403434 
    … …  
    405436      zfact1 = 0.5 * 0.25 
    406437 
    407       ! Save ua and va trends 
    408       IF( l_trddyn )   THEN 
    409          ztdua(:,:,:) = ua(:,:,:)  
    410          ztdva(:,:,:) = va(:,:,:)  
    411          zcu(:,:,:) = 0.e0 
    412          zcv(:,:,:) = 0.e0 
    413       ENDIF 
    414  
     438!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
     439!$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
    415440      !                                                ! =============== 
    416441      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    419444         ! Potential vorticity and horizontal fluxes 
    420445         ! ----------------------------------------- 
    421          IF( lk_sco ) THEN 
     446         SELECT CASE( cd_vor )      ! vorticity considered 
     447         CASE ( 'COR' )                   ! planetary vorticcity (Coriolis) 
     448            zwz(:,:) = ff(:,:) 
     449         CASE ( 'VOR' )                   ! relative  vorticity 
     450            zwz(:,:) = rotn(:,:,jk) 
     451         CASE ( 'TOT' )                   ! total     vorticity (planetary + relative) 
     452            zwz(:,:) = rotn(:,:,jk) + ff(:,:) 
     453         END SELECT 
     454 
     455         IF( ln_sco ) THEN 
    422456            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop  
    423457               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    424                   zwz(ji,jj,jk) = ( rotn(ji,jj,jk) + ff(ji,jj) ) / fse3f(ji,jj,jk) 
    425                   zwx(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 
    426                   zwy(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 
     458                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 
     459                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 
     460                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 
    427461               END DO 
    428462            END DO 
    … …  
    430464            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop  
    431465               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    432                   zwz(ji,jj,jk) = rotn(ji,jj,jk) + ff(ji,jj) 
    433                   zwx(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk) 
    434                   zwy(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk) 
     466                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 
     467                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 
    435468               END DO 
    436469            END DO 
    … …  
    442475         DO jj = 2, jpjm1 
    443476            DO ji = fs_2, fs_jpim1   ! vector opt. 
    444                zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1,jk) + zwy(ji+1,jj-1,jk)   & 
    445                                             + zwy(ji  ,jj  ,jk) + zwy(ji+1,jj  ,jk) ) 
    446                zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ,jk) + zwx(ji-1,jj+1,jk)   & 
    447                                             + zwx(ji  ,jj  ,jk) + zwx(ji  ,jj+1,jk) ) 
    448  
    449                zua  = zuav * ( zwz(ji  ,jj-1,jk) + zwz(ji,jj,jk) ) 
    450                zva  = zvau * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj,jk) ) 
    451  
    452                ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    453                va(ji,jj,jk) = va(ji,jj,jk) + zva 
     477               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   & 
     478                                            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) 
     479               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   & 
     480                                            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) 
     481 
     482               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     483               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    454484            END DO   
    455485         END DO   
    … …  
    458488      !                                                ! =============== 
    459489 
    460  
    461       ! save the relative & planetary vorticity trends for diagnostic 
    462       ! momentum trends 
    463       IF( l_trddyn )   THEN 
    464          ! Compute the planetary vorticity term trend 
    465          !                                                ! =============== 
    466          DO jk = 1, jpkm1                                 ! Horizontal slab 
    467             !                                             ! =============== 
    468             DO jj = 2, jpjm1 
    469                DO ji = fs_2, fs_jpim1   ! vector opt. 
    470                   zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1,jk) + zwy(ji+1,jj-1,jk)   & 
    471                      &                         + zwy(ji  ,jj  ,jk) + zwy(ji+1,jj  ,jk) ) 
    472                   zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ,jk) + zwx(ji-1,jj+1,jk)   & 
    473                      &                         + zwx(ji  ,jj  ,jk) + zwx(ji  ,jj+1,jk) ) 
    474 # if defined key_s_coord 
    475                   zcu(ji,jj,jk)  = zuav * ( ff(ji  ,jj-1) / fse3f(ji  ,jj-1,jk)   & 
    476                     &                     + ff(ji  ,jj  ) / fse3f(ji  ,jj  ,jk) ) 
    477                   zcv(ji,jj,jk)  = zvau * ( ff(ji-1,jj  ) / fse3f(ji-1,jj  ,jk)   & 
    478                     &                     + ff(ji  ,jj  ) / fse3f(ji  ,jj  ,jk) ) 
    479 # else 
    480                   zcu(ji,jj,jk) = zuav * ( ff(ji  ,jj-1) + ff(ji,jj) ) 
    481                   zcv(ji,jj,jk) = zvau * ( ff(ji-1,jj  ) + ff(ji,jj) ) 
    482 # endif 
    483                END DO   
    484             END DO   
    485             !                                             ! =============== 
    486          END DO                                           !   End of slab 
    487          !                                                ! =============== 
    488  
    489          ! Compute the relative vorticity term trend 
    490          ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:)  
    491          ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:)  
    492  
    493          CALL trd_mod(zcu  , zcv  , jpdtdpvo, 'DYN', kt) 
    494          CALL trd_mod(zcu  , zcv  , jpdtddat, 'DYN', kt) 
    495          CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 
    496       ENDIF 
    497  
    498       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    499          CALL prt_ctl(tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, & 
    500             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    501       ENDIF 
    502  
    503    END SUBROUTINE dyn_vor_enstrophy 
    504  
    505  
    506    SUBROUTINE dyn_vor_ene_ens( kt ) 
    507       !!---------------------------------------------------------------------- 
    508       !!                ***  ROUTINE dyn_vor_ene_ens  *** 
     490   END SUBROUTINE vor_ens 
     491 
     492 
     493   SUBROUTINE vor_een( kt, cd_vor, pua, pva ) 
     494      !!---------------------------------------------------------------------- 
     495      !!                ***  ROUTINE vor_een  *** 
    509496      !! 
    510497      !! ** Purpose :   Compute the now total vorticity trend and add it to  
    … …  
    516503      !!      when horizontal divergence is zero. 
    517504      !!      The trend of the vorticity term is given by: 
    518       !!       * s-coordinate (lk_sco=T), the e3. are inside the derivatives: 
     505      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 
    519506      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 
    520507      !!      Add this trend to the general momentum trend (ua,va): 
    … …  
    534521      !!   9.0  !  04-08  (C. Talandier)  New trends organization 
    535522      !!---------------------------------------------------------------------- 
    536       !! * Modules used      
    537       USE oce, ONLY :    ztdua => ta,    & ! use ta as 3D workspace    
    538                          ztdva => sa       ! use sa as 3D workspace    
    539  
    540523      !! * Arguments 
    541524      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index 
     525      CHARACTER(len=3) , INTENT( in ) ::   & 
     526         cd_vor        ! define the vorticity considered 
     527         !             ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 
     528      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   & 
     529         pua, pva                      ! ???? 
    542530 
    543531      !! * Local declarations 
    … …  
    547535      REAL(wp), DIMENSION(jpi,jpj) ::   & 
    548536         zwx, zwy, zwz,                 &  ! temporary workspace 
    549          ztnw, ztne, ztsw, ztse,        &  !    "           " 
    550          zcor                              ! potential planetary vorticity (f/e3) 
    551       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    552          zcu, zcv                          ! temporary workspace   
     537         ztnw, ztne, ztsw, ztse            !    "           " 
    553538      REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE ::   & 
    554539         ze3f 
    … …  
    557542      IF( kt == nit000 ) THEN 
    558543         IF(lwp) WRITE(numout,*) 
    559          IF(lwp) WRITE(numout,*) 'dyn_vor_ene_ens : vorticity term: energy and enstrophy conserving scheme' 
    560          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 
     544         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 
     545         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    561546 
    562547         DO jk = 1, jpk 
    … …  
    565550                  ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    566551                     &             + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * 0.25 
    567 !!!               ze3f(ji,jj,jk) = MAX( ze3f(ji,jj,jk) , 1.e-20) 
    568552                  IF( ze3f(ji,jj,jk) /= 0.e0 )   ze3f(ji,jj,jk) = 1.e0 / ze3f(ji,jj,jk) 
    569553               END DO 
    … …  
    576560      zfac12 = 1.e0 / 12.e0 
    577561 
    578       ! Save ua and va trends 
    579       IF( l_trddyn )   THEN 
    580          ztdua(:,:,:) = ua(:,:,:)  
    581          ztdva(:,:,:) = va(:,:,:)  
    582          zcu(:,:,:) = 0.e0 
    583          zcv(:,:,:) = 0.e0 
    584       ENDIF 
    585562       
     563!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 
     564!$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 
    586565      !                                                ! =============== 
    587566      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    590569         ! Potential vorticity and horizontal fluxes 
    591570         ! ----------------------------------------- 
    592 !!!bug   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) / fse3f(:,:,jk) 
    593          zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 
     571         SELECT CASE( cd_vor )      ! vorticity considered 
     572         CASE ( 'COR' )                   ! planetary vorticcity (Coriolis) 
     573            zwz(:,:) = ff(:,:) * ze3f(:,:,jk) 
     574         CASE ( 'VOR' )                   ! relative  vorticity 
     575            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 
     576         CASE ( 'TOT' )                   ! total     vorticity (planetary + relative) 
     577            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 
     578         END SELECT 
     579 
    594580         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    595581         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    596          zcor(:,:) = ff(:,:) * ze3f(:,:,jk) 
    597582 
    598583         ! Compute and add the vorticity term trend 
    … …  
    621606               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    622607                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    623                ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    624                va(ji,jj,jk) = va(ji,jj,jk) + zva 
     608               pua(ji,jj,jk) = pua(ji,jj,jk) + zua 
     609               pva(ji,jj,jk) = pva(ji,jj,jk) + zva 
    625610            END DO   
    626611         END DO   
    … …  
    629614      !                                                ! =============== 
    630615 
    631       ! save the relative & planetary vorticity trends for diagnostic 
    632       ! momentum trends 
    633       IF( l_trddyn )   THEN 
    634          ! Compute the planetary vorticity term trend 
    635          !                                                ! =============== 
    636          DO jk = 1, jpkm1                                 ! Horizontal slab 
    637             !                                             ! =============== 
    638             DO jj = 2, jpjm1 
    639                DO ji = fs_2, fs_jpim1   ! vector opt. 
    640                   zcu(ji,jj,jk) = + zfac12 / e1u(ji,jj) * (  zcor(ji,jj  ) * zwy(ji  ,jj  ) + zcor(ji+1,jj) * zwy(ji+1,jj  )   & 
    641                     &                                      + zcor(ji,jj  ) * zwy(ji  ,jj-1) + zcor(ji+1,jj) * zwy(ji+1,jj-1) ) 
    642                   zcv(ji,jj,jk) = - zfac12 / e2v(ji,jj) * (  zcor(ji,jj+1) * zwx(ji-1,jj+1) + zcor(ji,jj+1) * zwx(ji  ,jj+1)   & 
    643                     &                                      + zcor(ji,jj  ) * zwx(ji-1,jj  ) + zcor(ji,jj  ) * zwx(ji  ,jj  ) ) 
    644                END DO   
    645             END DO   
    646             !                                             ! =============== 
    647          END DO                                           !   End of slab 
    648          !                                                ! =============== 
    649  
    650          ! Compute the relative vorticity term trend 
    651          ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:)  
    652          ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:)  
    653  
    654          CALL trd_mod(zcu  , zcv  , jpdtdpvo, 'DYN', kt) 
    655          CALL trd_mod(zcu  , zcv  , jpdtddat, 'DYN', kt) 
    656          CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 
    657       ENDIF 
    658  
    659       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    660          CALL prt_ctl(tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, & 
    661             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    662       ENDIF 
    663  
    664    END SUBROUTINE dyn_vor_ene_ens 
    665  
    666  
    667    SUBROUTINE dyn_vor_ctl 
     616   END SUBROUTINE vor_een 
     617 
     618 
     619   SUBROUTINE vor_ctl 
    668620      !!--------------------------------------------------------------------- 
    669       !!                  ***  ROUTINE dyn_vor_ctl  *** 
     621      !!                  ***  ROUTINE vor_ctl  *** 
    670622      !! 
    671623      !! ** Purpose :   Control the consistency between cpp options for 
    … …  
    691643      IF(lwp) THEN 
    692644         WRITE(numout,*) 
    693          WRITE(numout,*) 'dyn_vor_ctl : vorticity term : read namelist and control the consistency' 
     645         WRITE(numout,*) 'dyn:vor_ctl : vorticity term : read namelist and control the consistency' 
    694646         WRITE(numout,*) '~~~~~~~~~~~' 
    695647         WRITE(numout,*) '              Namelist nam_dynvor : oice of the vorticity term scheme' 
     648         WRITE(numout,*) '                 energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene 
    696649         WRITE(numout,*) '                 enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens 
    697          WRITE(numout,*) '                 energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene 
    698650         WRITE(numout,*) '                 mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix 
    699651         WRITE(numout,*) '                 enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een 
    … …  
    701653 
    702654      ioptio = 0 
    703  
     655      IF( ln_dynvor_ene ) THEN 
     656         nvor = 0 
     657         IF(lwp) WRITE(numout,*) 
     658         IF(lwp) WRITE(numout,*) '              vorticity term : energy conserving scheme' 
     659         ioptio = ioptio + 1 
     660      ENDIF 
    704661      IF( ln_dynvor_ens ) THEN 
     662         nvor = 1 
    705663         IF(lwp) WRITE(numout,*) 
    706664         IF(lwp) WRITE(numout,*) '              vorticity term : enstrophy conserving scheme' 
    707665         ioptio = ioptio + 1 
    708666      ENDIF 
    709       IF( ln_dynvor_ene ) THEN 
    710          IF(lwp) WRITE(numout,*) 
    711          IF(lwp) WRITE(numout,*) '              vorticity term : energy conserving scheme' 
    712          ioptio = ioptio + 1 
    713       ENDIF 
    714667      IF( ln_dynvor_mix ) THEN 
     668         nvor = 2 
    715669         IF(lwp) WRITE(numout,*) 
    716670         IF(lwp) WRITE(numout,*) '              vorticity term : mixed enstrophy/energy conserving scheme' 
    … …  
    718672      ENDIF 
    719673      IF( ln_dynvor_een ) THEN 
     674         nvor = 3 
    720675         IF(lwp) WRITE(numout,*) 
    721676         IF(lwp) WRITE(numout,*) '              vorticity term : energy and enstrophy conserving scheme' 
    722677         ioptio = ioptio + 1 
    723678      ENDIF 
    724       IF ( ioptio /= 1 .AND. .NOT. lk_esopa ) THEN 
    725           WRITE(numout,cform_err) 
    726           IF(lwp) WRITE(numout,*) ' use ONE and ONLY one vorticity scheme' 
    727           nstop = nstop + 1 
    728       ENDIF 
    729  
    730    END SUBROUTINE dyn_vor_ctl 
     679      IF( ioptio /= 1 .AND. .NOT. lk_esopa ) THEN 
     680         WRITE(numout,cform_err) 
     681         IF(lwp) WRITE(numout,*) ' use ONE and ONLY one vorticity scheme' 
     682         nstop = nstop + 1 
     683      ENDIF 
     684      IF( lk_esopa ) THEN 
     685         nvor = -1 
     686         IF(lwp ) WRITE(numout,*) '          esopa test: use all lateral physics options' 
     687      ENDIF 
     688      IF(lwp) WRITE(numout,*) '             choice of vor_...                nvor   = ', nvor 
     689 
     690   END SUBROUTINE vor_ctl 
    731691 
    732692!!============================================================================== 

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

    r258 r455  
    3434CONTAINS 
    3535 
    36 #if defined key_autotasking 
    37    !!---------------------------------------------------------------------- 
    38    !!   'key_autotasking'                              j-k-i loops (j-slab) 
     36#if defined key_mpp_omp 
     37   !!---------------------------------------------------------------------- 
     38   !!   'key_mpp_omp'        OpenMP / NEC autotasking: j-k-i loops (j-slab) 
    3939   !!---------------------------------------------------------------------- 
    4040 
    … …  
    4646      !!      add it to the general trend of momentum equation. 
    4747      !! 
    48       !! ** Method  :   Use j-slab (j-k-i loops) for auto-tasking 
     48      !! ** Method  :   Use j-slab (j-k-i loops) for OpenMP / NEC autotasking 
    4949      !!      The now vertical advection of momentum is given by: 
    5050      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ] 

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

    r258 r455  
    1616   USE in_out_manager  ! I/O manager 
    1717   USE taumod          ! surface ocean stress 
    18    USE trdmod          ! ocean dynamics trends  
    19    USE trdmod_oce      ! ocean variables trends 
    2018   USE prtctl          ! Print control 
    2119 
    … …  
    6866      !! * Local declarations 
    6967      INTEGER ::   & 
    70          ji, jj, jk, jl,                 & ! dummy loop indices 
    71          ikbu, ikbum1 , ikbv, ikbvm1       ! temporary integers 
     68         ji, jj, jk, jl                    ! dummy loop indices 
    7269      REAL(wp) ::   & 
    7370         zrau0r, zlavmr, z2dt, zua, zva    ! temporary scalars 
    7471      REAL(wp), DIMENSION(jpi,jpk) ::    & 
    7572         zwx, zwy, zwz, zww                ! temporary workspace arrays 
    76       REAL(wp), DIMENSION(jpi,jpj) ::    & 
    77          ztsx, ztsy, ztbx, ztby            ! temporary workspace arrays 
    7873      !!---------------------------------------------------------------------- 
    7974 
    … …  
    8984      zlavmr = 1. / float( n_zdfexp )                      ! inverse of the number of sub time step 
    9085      z2dt = 2. * rdt                                      ! Leap-frog environnement 
    91       ztsx(:,:) = 0.e0 
    92       ztsy(:,:) = 0.e0  
    93       ztbx(:,:) = 0.e0 
    94       ztby(:,:) = 0.e0 
    95  
    96       ! Save ua and va trends 
    97       IF( l_trddyn )   THEN 
    98          ztdua(:,:,:) = ua(:,:,:)  
    99          ztdva(:,:,:) = va(:,:,:)  
    100       ENDIF 
    10186 
    10287      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt    ! Euler time stepping when starting from rest 
    … …  
    150135      !                                                ! =============== 
    151136 
    152       ! save the vertical diffusion trends for diagnostic 
    153       ! momentum trends 
    154       IF( l_trddyn )  THEN  
    155          ! save the total vertical momentum diffusive trend 
    156          ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 
    157          ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) 
    158   
    159          ! subtract and save surface and momentum fluxes 
    160          !                                                ! =============== 
    161          DO jj = 2, jpjm1                                 !  Horizontal slab 
    162             !                                             ! =============== 
    163             DO ji = 2, jpim1 
    164                ! save the surface momentum fluxes  
    165                ztsx(ji,jj) = zwy(ji,1) / fse3u(ji,jj,1) 
    166                ztsy(ji,jj) = zww(ji,1) / fse3v(ji,jj,1) 
    167                ! save bottom friction momentum fluxes  
    168                ikbu   = MIN( mbathy(ji+1,jj), mbathy(ji,jj) ) 
    169                ikbum1 = MAX( ikbu-1, 1 ) 
    170                ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    171                ikbvm1 = MAX( ikbv-1, 1 ) 
    172                ztbx(ji,jj) = avmu(ji,jj,ikbu) * zwx(ji,ikbum1)   & 
    173                                / ( fse3u(ji,jj,ikbum1) * fse3uw(ji,jj,ikbu) ) 
    174                ztby(ji,jj) = avmv(ji,jj,ikbv) * zwz(ji,ikbvm1)   & 
    175                                / ( fse3v(ji,jj,ikbvm1) * fse3vw(ji,jj,ikbv) ) 
    176                ! subtract surface forcing and bottom friction trend from vertical 
    177                ! diffusive momentum trend 
    178                ztdua(ji,jj,1     ) = ztdua(ji,jj,1     ) - ztsx(ji,jj) 
    179                ztdua(ji,jj,ikbum1) = ztdua(ji,jj,ikbum1) - ztbx(ji,jj) 
    180                ztdva(ji,jj,1     ) = ztdva(ji,jj,1     ) - ztsy(ji,jj) 
    181                ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj) 
    182             END DO 
    183             !                                             ! =============== 
    184          END DO                                           !   End of slab 
    185          !                                                ! =============== 
    186  
    187          CALL trd_mod(ztdua, ztdva, jpdtdzdf, 'DYN', kt) 
    188          ztdua(:,:,:) = 0.e0 
    189          ztdva(:,:,:) = 0.e0 
    190          ztdua(:,:,1) = ztsx(:,:) 
    191          ztdva(:,:,1) = ztsy(:,:) 
    192          CALL trd_mod(ztdua , ztdva , jpdtdswf, 'DYN', kt) 
    193          ztdua(:,:,:) = 0.e0 
    194          ztdva(:,:,:) = 0.e0 
    195          ztdua(:,:,1) = ztbx(:,:) 
    196          ztdva(:,:,1) = ztby(:,:) 
    197          CALL trd_mod(ztdua , ztdva , jpdtdbfr, 'DYN', kt) 
    198       ENDIF 
    199  
    200       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    201          CALL prt_ctl(tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask, & 
    202             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    203       ENDIF 
    204  
    205137   END SUBROUTINE dyn_zdf_exp 
    206138 

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

    r258 r455  
    2020   USE in_out_manager  ! I/O manager 
    2121   USE taumod          ! surface ocean stress 
    22    USE trdmod          ! ocean dynamics trends  
    23    USE trdmod_oce      ! ocean variables trends 
    2422   USE prtctl          ! Print control 
    2523 
    … …  
    6058      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive 
    6159      !!               mixing trend. 
    62       !!             - Save the trends in (ztdua,ztdva) ('l_trddyn') 
    6360      !! 
    6461      !! History : 
    … …  
    7673 
    7774      !! * Local declarations 
    78       INTEGER ::   & 
    79          ji, jj, jk,                  &  ! dummy loop indices 
    80          ikbu, ikbum1, ikbv, ikbvm1      ! temporary integers 
     75      INTEGER ::   ji, jj, jk            ! dummy loop indices 
    8176      REAL(wp) ::   & 
    8277         zrau0r, z2dt,                &  ! temporary scalars 
    8378         z2dtf, zcoef, zzws, zrhs        !    "         " 
    84       REAL(wp), DIMENSION(jpi,jpj) ::   & 
    85          ztsx, ztsy, ztbx, ztby          ! temporary workspace arrays 
    8679      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    87          zwi, ztdua, ztdva               ! temporary workspace arrays 
     80         zwi                             ! temporary workspace arrays 
    8881      !!---------------------------------------------------------------------- 
    8982 
    … …  
    9891      zrau0r = 1. / rau0      ! inverse of the reference density 
    9992      z2dt   = 2. * rdt       ! Leap-frog environnement 
    100       ztsx(:,:) = 0.e0 
    101       ztsy(:,:) = 0.e0  
    102       ztbx(:,:) = 0.e0 
    103       ztby(:,:) = 0.e0 
     93 
    10494      ! Euler time stepping when starting from rest 
    10595      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
    106  
    107       ! Save previous ua and va trends 
    108       IF( l_trddyn )   THEN 
    109          ztdua(:,:,:) = ua(:,:,:)  
    110          ztdva(:,:,:) = va(:,:,:)  
    111       ENDIF 
    11296 
    11397      ! 1. Vertical diffusion on u 
    … …  
    194178         END DO 
    195179      END DO 
    196  
    197       IF( l_trddyn )  THEN  
    198          ! diagnose surface and bottom momentum fluxes 
    199          DO jj = 2, jpjm1    
    200             DO ji = fs_2, fs_jpim1   ! vector opt. 
    201                ! save the surface forcing momentum fluxes 
    202                ztsx(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) 
    203                ! save bottom friction momentum fluxes 
    204                ikbu   = MIN( mbathy(ji+1,jj), mbathy(ji,jj) ) 
    205                ikbum1 = MAX( ikbu-1, 1 ) 
    206                ztbx(ji,jj) = - avmu(ji,jj,ikbu) * ua(ji,jj,ikbum1)   & 
    207                   / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) ) 
    208                ! subtract surface forcing and bottom friction trend from vertical 
    209                ! diffusive momentum trend 
    210                ztdua(ji,jj,1     ) = ztdua(ji,jj,1     ) - ztsx(ji,jj) 
    211                ztdua(ji,jj,ikbum1) = ztdua(ji,jj,ikbum1) - ztbx(ji,jj) 
    212             END DO 
    213          END DO 
    214       ENDIF 
    215180 
    216181      ! Normalization to obtain the general momentum trend ua 
    … …  
    309274      END DO 
    310275 
    311       IF( l_trddyn )  THEN  
    312          ! diagnose surface and bottom momentum fluxes 
    313          DO jj = 2, jpjm1    
    314             DO ji = fs_2, fs_jpim1   ! vector opt. 
    315                ! save the surface forcing momentum fluxes 
    316                ztsy(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) 
    317                ! save bottom friction momentum fluxes 
    318                ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    319                ikbvm1 = MAX( ikbv-1, 1 ) 
    320                ztby(ji,jj) = - avmv(ji,jj,ikbv) * va(ji,jj,ikbvm1)   & 
    321                   / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) ) 
    322                ! subtract surface forcing and bottom friction trend from vertical 
    323                ! diffusive momentum trend 
    324                ztdva(ji,jj,1     ) = ztdva(ji,jj,1     ) - ztsy(ji,jj) 
    325                ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj) 
    326             END DO 
    327          END DO 
    328       ENDIF 
     276! flux de surface doit etre calcule dans trdmod et boootom stress 
     277! deduit par integration verticale dans trmod pour jpdtdzdf 
     278!RB      IF( l_trddyn )  THEN  
     279!         ! diagnose surface and bottom momentum fluxes 
     280!         DO jj = 2, jpjm1    
     281!            DO ji = fs_2, fs_jpim1   ! vector opt. 
     282!               ! save the surface forcing momentum fluxes 
     283!               ztsy(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) 
     284!               ztsx(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) 
     285!               ! save bottom friction momentum fluxes 
     286!               ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
     287!               ikbvm1 = MAX( ikbv-1, 1 ) 
     288!               ztby(ji,jj) = - avmv(ji,jj,ikbv) * va(ji,jj,ikbvm1)   & 
     289!                  / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) ) 
     290!               ! subtract surface forcing and bottom friction trend from vertical 
     291!               ! diffusive momentum trend 
     292!               ztdva(ji,jj,1     ) = ztdva(ji,jj,1     ) - ztsy(ji,jj) 
     293!               ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj) 
     294!            END DO 
     295!         END DO 
     296!      ENDIF 
    329297 
    330298      ! Normalization to obtain the general momentum trend va 
    … …  
    337305      END DO 
    338306 
    339       ! save the vertical diffusion trends for diagnostic 
    340       ! momentum trends 
    341       IF( l_trddyn )  THEN  
    342          ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) 
    343          ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) 
    344  
    345          CALL trd_mod(ztdua, ztdva, jpdtdzdf, 'DYN', kt) 
    346          ztdua(:,:,:) = 0.e0 
    347          ztdva(:,:,:) = 0.e0 
    348          ztdua(:,:,1) = ztsx(:,:) 
    349          ztdva(:,:,1) = ztsy(:,:) 
    350          CALL trd_mod(ztdua , ztdva , jpdtdswf, 'DYN', kt) 
    351          ztdua(:,:,:) = 0.e0 
    352          ztdva(:,:,:) = 0.e0 
    353          ztdua(:,:,1) = ztbx(:,:) 
    354          ztdva(:,:,1) = ztby(:,:) 
    355          CALL trd_mod(ztdua , ztdva , jpdtdbfr, 'DYN', kt) 
    356       ENDIF 
    357  
    358       IF(ln_ctl) THEN         ! print sum trends (used for debugging) 
    359          CALL prt_ctl(tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask, & 
    360             &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 
    361       ENDIF 
    362  
    363307   END SUBROUTINE dyn_zdf_imp 
    364308 

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

    r258 r455  
    2626CONTAINS 
    2727 
    28 #if defined key_autotasking 
     28#if defined key_mpp_omp 
    2929   !!---------------------------------------------------------------------- 
    30    !!   'key_autotasking'                               j-k-i loop (j-slab) 
     30   !!   'key_mpp_omp'                                   j-k-i loop (j-slab) 
    3131   !!---------------------------------------------------------------------- 
    3232 
    … …  
    6464         IF(lwp) WRITE(numout,*) 
    6565         IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.' 
    66          IF(lwp) WRITE(numout,*) '~~~~~~~   auto-tasking case : j-k-i loop ' 
     66         IF(lwp) WRITE(numout,*) '~~~~~~~                     j-k-i loops' 
    6767 
    6868         ! bottom boundary condition: w=0 (set once for all) 
Note: See TracChangeset for help on using the changeset viewer.