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

Changeset 9528


Ignore:
Timestamp:
2018-04-30T15:39:18+02:00 (4 years ago)
Author:
gm
Message:

dev_merge_2017: add to new vorticity scheme (EET= EEN using e3t instead of e3f, and ENT= energy conserving at T-point). NB: not used in ref config.

Location:
branches/2017/dev_merge_2017/NEMOGCM
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg

    r9527 r9528  
    167167&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection) 
    168168!----------------------------------------------------------------------- 
    169    ln_dynvor_ene = .true.  !  enstrophy conserving scheme 
     169   ln_dynvor_ene = .true.  !  energy conserving scheme 
    170170/ 
    171171!----------------------------------------------------------------------- 
  • branches/2017/dev_merge_2017/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg

    r9527 r9528  
    168168&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection) 
    169169!----------------------------------------------------------------------- 
    170    ln_dynvor_ene = .true.  !  enstrophy conserving scheme 
     170   ln_dynvor_ene = .true.  !  energy conserving scheme 
    171171/ 
    172172!----------------------------------------------------------------------- 
  • branches/2017/dev_merge_2017/NEMOGCM/CONFIG/SHARED/namelist_ref

    r9527 r9528  
    22!! NEMO/OCE :   Reference namelist_ref                                !! 
    33!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    4 !! NEMO/OPA  :  1 - run manager      (namrun) 
    5 !! namelists    2 - Domain           (namcfg, namdom, namtsd, namcrs, namc1d, namc1d_uvd) 
    6 !!              3 - Surface boundary (namsbc, namsbc_flx, namsbc_blk, namsbc_cpl, 
     4!! NEMO/OPA  :  1 - Domain & run manager (namrun, namcfg, namdom, namtsd, namcrs, namc1d, namc1d_uvd) 
     5!! namelists    2 - Surface boundary (namsbc, namsbc_flx, namsbc_blk, namsbc_cpl, 
    76!!                                    namsbc_sas, namtra_qsr, namsbc_rnf, 
    87!!                                    namsbc_isf, namsbc_iscpl, namsbc_apr,  
    98!!                                    namsbc_ssr, namsbc_wave, namberg) 
    10 !!              4 - lateral boundary (namlbc, namagrif, nambdy, nambdy_tide) 
    11 !!              5 - top/bot boundary (namdrg, namdrg_top, namdrg_bot, nambbc, nambbl) 
    12 !!              6 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_eiv, namtra_dmp) 
    13 !!              7 - dynamics         (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf) 
    14 !!              8 - Vertical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_iwm) 
    15 !!              9 - miscellaneous    (nammpp, namctl) 
    16 !!             10 - diagnostics      (namnc4, namtrd, namspr, namflo, namhsb, namsto) 
    17 !!             11 - Obs & Assim      (namobs, nam_asminc) 
     9!!              3 - lateral boundary (namlbc, namagrif, nambdy, nambdy_tide) 
     10!!              4 - top/bot boundary (namdrg, namdrg_top, namdrg_bot, nambbc, nambbl) 
     11!!              5 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_eiv, namtra_dmp) 
     12!!              6 - dynamics         (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf) 
     13!!              7 - Vertical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_iwm) 
     14!!              8 - diagnostics      (namnc4, namtrd, namspr, namflo, namhsb) 
     15!!              9 - Obs & Assim      (namobs, nam_asminc) 
     16!!             10 - miscellaneous    (nammpp, namctl, namsto) 
    1817!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1918 
     
    651650&namdrg        !   top/bottom drag coefficient                          (default: NO selection) 
    652651!----------------------------------------------------------------------- 
    653    ln_OFF     = .false.    !  free-slip       : Cd = 0                  (F => fill namdrg_bot 
    654    ln_lin     = .false.    !      linear  drag: Cd = Cd0 Uc0                   &   namdrg_top) 
    655    ln_non_lin = .false.    !  non-linear  drag: Cd = Cd0 |U| 
    656    ln_loglayer= .false.    !  logarithmic drag: Cd = vkarmn/log(z/z0) |U| 
    657    ! 
    658    ln_drgimp  = .true.     !  implicit top/bottom friction flag 
     652   ln_OFF      = .false.   !  free-slip       : Cd = 0                  (F => fill namdrg_bot 
     653   ln_lin      = .false.   !      linear  drag: Cd = Cd0 Uc0                   &   namdrg_top) 
     654   ln_non_lin  = .false.   !  non-linear  drag: Cd = Cd0 |U| 
     655   ln_loglayer = .false.   !  logarithmic drag: Cd = vkarmn/log(z/z0) |U| 
     656   ! 
     657   ln_drgimp   = .true.    !  implicit top/bottom friction flag 
    659658/ 
    660659!----------------------------------------------------------------------- 
    661660&namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T) 
    662661!----------------------------------------------------------------------- 
    663    rn_Cd0     =  1.e-3     !  drag coefficient [-] 
    664    rn_Uc0     =  0.4       !  ref. velocity [m/s] (linear drag=Cd0*Uc0)  
    665    rn_Cdmax   =  0.1       !  drag value maximum [-] (logarithmic drag) 
    666    rn_ke0     =  2.5e-3    !  background kinetic energy  [m2/s2] (non-linear cases) 
    667    rn_z0      =  3.0e-3    !  roughness [m] (ln_loglayer=T) 
    668    ln_boost   = .false.    !  =T regional boost of Cd0 ; =F constant 
    669      rn_boost =  50.          !  local boost factor  [-] 
     662   rn_Cd0      =  1.e-3    !  drag coefficient [-] 
     663   rn_Uc0      =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0)  
     664   rn_Cdmax    =  0.1      !  drag value maximum [-] (logarithmic drag) 
     665   rn_ke0      =  2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases) 
     666   rn_z0       =  3.0e-3   !  roughness [m] (ln_loglayer=T) 
     667   ln_boost    = .false.   !  =T regional boost of Cd0 ; =F constant 
     668      rn_boost =  50.         !  local boost factor  [-] 
    670669/ 
    671670!----------------------------------------------------------------------- 
    672671&namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
    673672!----------------------------------------------------------------------- 
    674    rn_Cd0     =  1.e-3    !  drag coefficient [-] 
    675    rn_Uc0     =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0)  
    676    rn_Cdmax   =  0.1      !  drag value maximum [-] (logarithmic drag) 
    677    rn_ke0     =  2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases) 
    678    rn_z0      =  3.e-3    !  roughness [m] (ln_loglayer=T) 
    679    ln_boost   = .false.   !  =T regional boost of Cd0 ; =F constant 
    680       rn_boost=  50.         !  local boost factor  [-] 
     673   rn_Cd0      =  1.e-3    !  drag coefficient [-] 
     674   rn_Uc0      =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0)  
     675   rn_Cdmax    =  0.1      !  drag value maximum [-] (logarithmic drag) 
     676   rn_ke0      =  2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases) 
     677   rn_z0       =  3.e-3    !  roughness [m] (ln_loglayer=T) 
     678   ln_boost    = .false.   !  =T regional boost of Cd0 ; =F constant 
     679      rn_boost =  50.         !  local boost factor  [-] 
    681680/ 
    682681!----------------------------------------------------------------------- 
     
    684683!----------------------------------------------------------------------- 
    685684   ln_trabbc   = .false.   !  Apply a geothermal heating at the ocean bottom 
    686      nn_geoflx     = 2       !  geothermal heat flux: = 1 constant flux 
    687       !                      !                        = 2 read variable flux [mW/m2] 
    688      rn_geoflx_cst = 86.4e-3 !  Constant value of geothermal heat flux       [mW/m2] 
     685      nn_geoflx     = 2       !  geothermal heat flux: = 1 constant flux 
     686      !                       !                        = 2 read variable flux [mW/m2] 
     687      rn_geoflx_cst = 86.4e-3 !  Constant value of geothermal heat flux       [mW/m2] 
    689688 
    690689   cn_dir      = './'      !  root directory for the geothermal data location 
     
    866865   ln_dynvor_ens = .false. !  enstrophy conserving scheme 
    867866   ln_dynvor_mix = .false. !  mixed scheme 
     867   ln_dynvor_enT = .false. !  energy conserving scheme (T-point) 
     868   ln_dynvor_eeT = .false. !  energy conserving scheme (een using e3t) 
    868869   ln_dynvor_een = .false. !  energy & enstrophy scheme 
    869       nn_een_e3f = 1          ! =0   e3f = mean masked e3t divided by 4 
    870       !                       ! =1   e3f = mean masked e3t divided by the sum of mask 
    871    ln_dynvor_msk = .false. !  vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)  ! PLEASE DO NOT ACTIVATE 
     870      nn_een_e3f = 1          ! =0  e3f = mi(mj(e3t))/4  
     871      !                       ! =1  e3f = mi(mj(e3t))/mi(mj( tmask)) 
     872   ln_dynvor_msk = .false. !  vorticity multiplied by fmask (=T)        ==>>> PLEASE DO NOT ACTIVATE 
     873      !                    !  (f-point vorticity schemes only) 
    872874/ 
    873875!----------------------------------------------------------------------- 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r9506 r9528  
    3636   USE sbcapr          ! surface boundary condition: atmospheric pressure 
    3737   USE dynadv    , ONLY: ln_dynadv_vec 
     38   USE dynvor          ! vortivity scheme indicators 
    3839   USE phycst          ! physical constants 
    3940   USE dynvor          ! vorticity term 
     
    103104      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 
    104105      ! 
    105       IF( ln_dynvor_een )   ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
    106          &                            ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 
     106      IF( ln_dynvor_een .OR. ln_dynvor_eeT )   & 
     107         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
     108         &               ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 
    107109         ! 
    108110      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) ) 
     
    149151      INTEGER  ::   ikbu, iktu, noffset   ! local integers 
    150152      INTEGER  ::   ikbv, iktv            !   -      - 
    151       REAL(wp) ::   r1_2dt_b, z2dt_bf        ! local scalars 
    152       REAL(wp) ::   zx1, zx2, zu_spg, zhura  !   -      - 
    153       REAL(wp) ::   zy1, zy2, zv_spg, zhvra  !   -      - 
    154       REAL(wp) ::   za0, za1, za2, za3       !   -      - 
    155       REAL(wp) ::   zmdi, zztmp              !   -      - 
     153      REAL(wp) ::   r1_2dt_b, z2dt_bf               ! local scalars 
     154      REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      - 
     155      REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      - 
     156      REAL(wp) ::   za0, za1, za2, za3              !   -      - 
     157      REAL(wp) ::   zmdi, zztmp            , z1_ht  !   -      - 
    156158      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 
    157159      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 
    158160      REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 
    159       REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e 
     161      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 
    160162      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 
    161163      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points 
     
    237239      ! 
    238240      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 
    239          IF( ln_dynvor_een ) THEN               !==  EEN scheme  ==! 
     241         ! 
     242         SELECT CASE( nvor_scheme ) 
     243         CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme) 
    240244            SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
    241245            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     
    250254               DO jj = 1, jpjm1 
    251255                  DO ji = 1, jpim1 
    252                      zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     & 
    253                         &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   & 
    254                         &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
    255                         &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
     256!!gm bug for ocean cavities 
     257!                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     & 
     258!                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   & 
     259!                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
     260!                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
     261!! replace by : 
     262                     zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      & 
     263                        &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   & 
     264                        &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      & 
     265                        &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   ) 
     266!!gm end 
    256267                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    257268                  END DO 
     
    270281            END DO 
    271282            ! 
    272          ELSE                                !== all other schemes (ENE, ENS, MIX) 
     283         CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme) 
     284            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     285            DO jj = 2, jpj 
     286               DO ji = 2, jpi 
     287                  z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 
     288                  ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht 
     289                  ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht 
     290                  ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 
     291                  ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht 
     292               END DO 
     293            END DO 
     294            ! 
     295         CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT ! 
     296            ! 
    273297            zwz(:,:) = 0._wp 
    274298            zhf(:,:) = 0._wp 
     
    280304!!gm  
    281305!!             
    282               IF( .NOT.ln_sco ) THEN 
     306            IF( .NOT.ln_sco ) THEN 
    283307   
    284308   !!gm  agree the JC comment  : this should be done in a much clear way 
     
    290314   !              ENDIF 
    291315   !              zhf(:,:) = gdepw_0(:,:,jk+1) 
    292                ELSE 
    293                  !zhf(:,:) = hbatf(:,:) 
    294                  DO jj = 1, jpjm1 
    295                    DO ji = 1, jpim1 
    296                      zhf(ji,jj) = MAX( 0._wp,                                & 
    297                                 & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     & 
    298                                 &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     & 
    299                                 &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     & 
    300                                 &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   & 
    301                                 & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +& 
    302                                 &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +& 
    303                                 &   rsmall  ) ) 
    304                    END DO 
    305                  END DO 
    306               END IF 
    307    
    308               DO jj = 1, jpjm1 
    309                  zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
    310               END DO 
     316               ! 
     317            ELSE 
     318               ! 
     319               !zhf(:,:) = hbatf(:,:) 
     320               DO jj = 1, jpjm1 
     321                  DO ji = 1, jpim1 
     322!!gm Bug here in case of ocean cavities, further more ht_0 is masked by definition ==>> remove mask multiplication  
     323!                     zhf(ji,jj) = MAX( 0._wp,                                & 
     324!                                & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     & 
     325!                                &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     & 
     326!                                &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     & 
     327!                                &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   & 
     328!                                & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +& 
     329!                                &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +& 
     330!                                &   rsmall  ) ) 
     331!!gm replaced by : 
     332                     zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          & 
     333                        &              + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      & 
     334                        &       / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          & 
     335                        &              + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  ) 
    311336!!gm end 
    312  
     337                  END DO 
     338               END DO 
     339            ENDIF 
     340            ! 
     341            DO jj = 1, jpjm1 
     342               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     343            END DO 
     344            ! 
    313345            DO jk = 1, jpkm1 
    314346               DO jj = 1, jpjm1 
     
    324356            END DO 
    325357            zwz(:,:) = ff_f(:,:) * zwz(:,:) 
    326          ENDIF 
     358         END SELECT 
    327359      ENDIF 
    328360      ! 
     
    369401      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    370402      !                                   ! -------------------------------------------------------- 
     403      ! 
    371404      zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes  
    372405      zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 
    373406      ! 
    374       IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
     407      SELECT CASE( nvor_scheme ) 
     408      CASE( np_ENT )                ! enstrophy conserving scheme (f-point) 
     409         DO jj = 2, jpjm1 
     410            DO ji = 2, jpim1   ! vector opt. 
     411               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    & 
     412                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   & 
     413                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   ) 
     414                  ! 
     415               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    & 
     416                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   &  
     417                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   )  
     418            END DO   
     419         END DO   
     420         !          
     421      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX 
    375422         DO jj = 2, jpjm1 
    376423            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    385432         END DO 
    386433         ! 
    387       ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme 
     434      CASE( np_ENS )                ! enstrophy conserving scheme (f-point) 
    388435         DO jj = 2, jpjm1 
    389436            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    397444         END DO 
    398445         ! 
    399       ELSEIF ( ln_dynvor_een ) THEN  ! enstrophy and energy conserving scheme 
     446      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)          
    400447         DO jj = 2, jpjm1 
    401448            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    411458         END DO 
    412459         ! 
    413       ENDIF  
     460      END SELECT 
    414461      ! 
    415462      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    416463      !                                   ! ---------------------------------------------------- 
    417464      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    418         IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters 
    419            DO jj = 2, jpjm1 
    420               DO ji = 2, jpim1  
    421                 ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
    422                      &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
    423                      &    MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     465         IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters 
     466            DO jj = 2, jpjm1 
     467               DO ji = 2, jpim1  
     468                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
     469                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     470                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     471                     &                                                         > rn_wdmin1 + rn_wdmin2 
     472                  ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
     473                     &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     474                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     475                  IF(ll_tmp1) THEN 
     476                     zcpx(ji,jj) = 1.0_wp 
     477                  ELSEIF(ll_tmp2) THEN 
     478                     ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     479                     zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
     480                                 &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     481                     zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
     482                  ELSE 
     483                     zcpx(ji,jj) = 0._wp 
     484                  ENDIF 
     485                  ! 
     486                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                & 
     487                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
     488                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
    424489                     &                                                       > rn_wdmin1 + rn_wdmin2 
    425                 ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
    426                      &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
    427                      &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    428     
    429                 IF(ll_tmp1) THEN 
    430                   zcpx(ji,jj) = 1.0_wp 
    431                 ELSE IF(ll_tmp2) THEN 
    432                   ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    433                   zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    434                               &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
    435                   zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
    436  
    437                 ELSE 
    438                   zcpx(ji,jj) = 0._wp 
    439                 END IF 
    440           
    441                 ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                & 
    442                      &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
    443                      &    MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
    444                      &                                                       > rn_wdmin1 + rn_wdmin2 
    445                 ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
    446                      &    MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
    447                      &    MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     490                  ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
     491                     &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     492                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    448493   
    449                 IF(ll_tmp1) THEN 
    450                   zcpy(ji,jj) = 1.0_wp 
    451                 ELSE IF(ll_tmp2) THEN 
    452                   ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    453                   zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    454                               &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
    455                   zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 
    456  
    457                 ELSE 
    458                   zcpy(ji,jj) = 0._wp 
    459                 END IF 
    460               END DO 
    461            END DO 
    462   
    463            DO jj = 2, jpjm1 
    464               DO ji = 2, jpim1 
    465                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    466                         &                        * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
    467                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    468                         &                        * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth 
    469  
    470               END DO 
    471            END DO 
    472  
     494                  IF(ll_tmp1) THEN 
     495                     zcpy(ji,jj) = 1.0_wp 
     496                  ELSE IF(ll_tmp2) THEN 
     497                     ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     498                     zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
     499                        &             / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     500                     zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  ) 
     501                  ELSE 
     502                     zcpy(ji,jj) = 0._wp 
     503                  ENDIF 
     504               END DO 
     505            END DO 
     506            ! 
     507            DO jj = 2, jpjm1 
     508               DO ji = 2, jpim1 
     509                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     510                     &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
     511                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     512                     &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth 
     513               END DO 
     514            END DO 
     515            ! 
    473516         ELSE 
    474517            ! 
     
    675718      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
    676719      vn_adv(:,:) = 0._wp 
    677       !                                             ! ==================== ! 
    678  
    679       IF (ln_wd_dl) THEN 
     720      ! 
     721      IF( ln_wd_dl ) THEN 
    680722         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)  
    681723         zvwdmask(:,:) = 0._wp  !  
    682          zuwdav2(:,:) = 0._wp  
    683          zvwdav2(:,:) = 0._wp    
     724         zuwdav2 (:,:) = 0._wp  
     725         zvwdav2 (:,:) = 0._wp    
    684726      END IF  
    685727 
    686  
     728      !                                             ! ==================== ! 
    687729      DO jn = 1, icycle                             !  sub-time-step loop  ! 
    688730         !                                          ! ==================== ! 
     
    715757            ! set wetting & drying mask at tracer points for this barotropic sub-step  
    716758            IF ( ln_wd_dl ) THEN  
    717  
     759               ! 
    718760               IF ( ln_wd_dl_rmp ) THEN  
    719761                  DO jj = 1, jpj                                  
     
    736778                        ELSE  
    737779                           ztwdmask(ji,jj) = 0._wp 
    738                         END IF 
     780                        ENDIF 
    739781                     END DO 
    740782                  END DO 
    741                END IF  
    742  
    743             END IF  
    744             
    745  
     783               ENDIF  
     784               ! 
     785            ENDIF  
     786            ! 
    746787            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
    747788               DO ji = 2, fs_jpim1   ! Vector opt. 
     
    756797            CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 
    757798            ! 
    758             zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    759             zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     799            zhup2_e(:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
     800            zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 
     801            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 
    760802         ELSE 
    761             zhup2_e (:,:) = hu_n(:,:) 
    762             zhvp2_e (:,:) = hv_n(:,:) 
     803            zhup2_e(:,:) = hu_n(:,:) 
     804            zhvp2_e(:,:) = hv_n(:,:) 
     805            zhtp2_e(:,:) = ht_n(:,:) 
    763806         ENDIF 
    764807         !                                                !* after ssh 
     
    795838         ENDIF 
    796839#endif 
    797          IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
     840         IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    798841 
    799842         IF ( ln_wd_dl ) THEN  
    800  
    801  
    802 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells  
    803  
     843            ! 
     844            ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells  
     845            ! 
    804846            DO jj = 1, jpjm1                                  
    805847               DO ji = 1, jpim1    
     
    821863               END DO 
    822864            END DO 
    823  
    824  
    825          END IF     
     865            ! 
     866         ENDIF     
    826867          
    827868         ! Sum over sub-time-steps to compute advective velocities 
     
    896937         ENDIF 
    897938         ! 
    898          zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
    899           &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
     939         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   & 
     940            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
     941          
    900942         IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters 
    901943           DO jj = 2, jpjm1 
     
    917959                ELSE 
    918960                  zcpx(ji,jj) = 0._wp 
    919                 END IF 
    920           
     961                ENDIF 
     962                ! 
    921963                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
    922964                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            & 
     
    929971                IF(ll_tmp1) THEN 
    930972                  zcpy(ji,jj) = 1.0_wp 
    931                 ELSE IF(ll_tmp2) THEN 
     973                ELSEIF(ll_tmp2) THEN 
    932974                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here 
    933975                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 
     
    935977                ELSE 
    936978                  zcpy(ji,jj) = 0._wp 
    937                 END IF 
     979                ENDIF 
    938980              END DO 
    939981           END DO 
    940          END IF 
     982         ENDIF 
    941983         ! 
    942984         ! Compute associated depths at U and V points: 
     
    955997               END DO 
    956998            END DO 
    957  
     999            ! 
    9581000         ENDIF 
    9591001         ! 
     
    9651007         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    9661008         ! 
    967          IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN     !==  energy conserving or mixed scheme  ==! 
     1009         SELECT CASE( nvor_scheme ) 
     1010         CASE( np_ENT )             ! energy conserving scheme (t-point) 
     1011         DO jj = 2, jpjm1 
     1012            DO ji = 2, jpim1   ! vector opt. 
     1013 
     1014               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 
     1015               z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
     1016             
     1017               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   & 
     1018                  &               * (  e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) )   & 
     1019                  &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   ) 
     1020                  ! 
     1021               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    & 
     1022                  &               * (  e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) )   &  
     1023                  &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   )  
     1024            END DO   
     1025         END DO   
     1026         !          
     1027         CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point) 
    9681028            DO jj = 2, jpjm1 
    9691029               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    9771037            END DO 
    9781038            ! 
    979          ELSEIF ( ln_dynvor_ens ) THEN                   !==  enstrophy conserving scheme  ==! 
     1039         CASE( np_ENS )             ! enstrophy conserving scheme (f-point) 
    9801040            DO jj = 2, jpjm1 
    9811041               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    9891049            END DO 
    9901050            ! 
    991          ELSEIF ( ln_dynvor_een ) THEN                   !==  energy and enstrophy conserving scheme  ==! 
     1051         CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f) 
    9921052            DO jj = 2, jpjm1 
    9931053               DO ji = fs_2, fs_jpim1   ! vector opt. 
    994                   zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
    995                      &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
    996                      &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1) &  
    997                      &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    998                   zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &  
    999                      &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
    1000                      &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &  
    1001                      &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     1054                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  & 
     1055                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  & 
     1056                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  &  
     1057                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  ) 
     1058                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  &  
     1059                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  & 
     1060                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  &  
     1061                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  ) 
    10021062               END DO 
    10031063            END DO 
    10041064            !  
    1005          ENDIF 
     1065         END SELECT 
    10061066         ! 
    10071067         ! Add tidal astronomical forcing if defined 
     
    13411401      !! ** Purpose : Read or write time-splitting arrays in restart file 
    13421402      !!---------------------------------------------------------------------- 
    1343       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    1344       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    1345       ! 
    1346       INTEGER ::   id1, id2   ! local integers 
     1403      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     1404      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    13471405      !!---------------------------------------------------------------------- 
    13481406      ! 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r9190 r9528  
    66   !!====================================================================== 
    77   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code 
    8    !!            5.0  ! 1991-11  (G. Madec) vor_ene, vor_mix: Original code 
     8   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code 
    99   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays 
    1010   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module 
     
    1919   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 
    2020   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis 
     21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet) 
     22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation 
    2123   !!---------------------------------------------------------------------- 
    2224 
     
    5052 
    5153   !                                   !!* Namelist namdyn_vor: vorticity term 
    52    LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme    (ENE) 
    53    LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme (ENS) 
    54    LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                (MIX) 
    55    LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme (EEN) 
     54   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS) 
     55   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE) 
     56   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT) 
     57   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET) 
     58   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN) 
    5659   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
     60   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX) 
    5761   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
    5862 
    59    INTEGER ::   nvor_scheme        ! choice of the type of advection scheme 
    60    !                               ! associated indices: 
     63   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme 
     64   !                                       ! associated indices: 
     65   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme 
    6166   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme 
    62    INTEGER, PUBLIC, PARAMETER ::   np_ENS = 2   ! ENS scheme 
    63    INTEGER, PUBLIC, PARAMETER ::   np_MIX = 3   ! MIX scheme 
     67   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity) 
     68   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t) 
    6469   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme 
     70   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme 
    6571 
    6672   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
    6773   !                               ! associated indices: 
    68    INTEGER, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
    69    INTEGER, PARAMETER ::   np_RVO = 2         ! relative vorticity 
    70    INTEGER, PARAMETER ::   np_MET = 3         ! metric term 
    71    INTEGER, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity) 
    72    INTEGER, PARAMETER ::   np_CME = 5         ! Coriolis + metric term 
     74   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
     75   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity 
     76   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term 
     77   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity) 
     78   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term 
     79 
     80   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation 
     81   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -  
     82   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation 
     83   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -  
    7384    
    7485   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
     
    109120         ztrdv(:,:,:) = va(:,:,:) 
    110121         SELECT CASE( nvor_scheme ) 
     122         CASE( np_ENS )           ;   CALL vor_ens( kt, ncor, un , vn , ua, va )   ! enstrophy conserving scheme 
     123            IF( ln_stcor )            CALL vor_ens( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    111124         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, ncor, un , vn , ua, va )   ! energy conserving scheme 
    112125            IF( ln_stcor )            CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    113          CASE( np_ENS )           ;   CALL vor_ens( kt, ncor, un , vn , ua, va )   ! enstrophy conserving scheme 
    114             IF( ln_stcor )            CALL vor_ens( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     126         CASE( np_ENT )           ;   CALL vor_enT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (T-pts) 
     127            IF( ln_stcor )            CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     128         CASE( np_EET )           ;   CALL vor_eeT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (een with e3t) 
     129            IF( ln_stcor )            CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    115130         CASE( np_EEN )           ;   CALL vor_een( kt, ncor, un , vn , ua, va )   ! energy & enstrophy scheme 
    116131            IF( ln_stcor )            CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     
    124139            ztrdv(:,:,:) = va(:,:,:) 
    125140            SELECT CASE( nvor_scheme ) 
     141            CASE( np_ENT )           ;   CALL vor_enT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts) 
     142            CASE( np_EET )           ;   CALL vor_eeT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (een with e3t) 
    126143            CASE( np_ENE )           ;   CALL vor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme 
    127144            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, nrvm, un , vn , ua, va )  ! enstrophy conserving scheme 
     
    138155         ! 
    139156         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==! 
     157         CASE( np_ENT )                        !* energy conserving scheme  (T-pts) 
     158                             CALL vor_enT( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
     159            IF( ln_stcor )   CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     160         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t) 
     161                             CALL vor_eeT( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
     162            IF( ln_stcor )   CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    140163         CASE( np_ENE )                        !* energy conserving scheme 
    141164                             CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
     
    164187 
    165188 
     189   SUBROUTINE vor_enT( kt, kvor, pu, pv, pu_rhs, pv_rhs ) 
     190      !!---------------------------------------------------------------------- 
     191      !!                  ***  ROUTINE vor_enT  *** 
     192      !! 
     193      !! ** Purpose :   Compute the now total vorticity trend and add it to  
     194      !!      the general trend of the momentum equation. 
     195      !! 
     196      !! ** Method  :   Trend evaluated using now fields (centered in time)  
     197      !!       and t-point evaluation of vorticity (planetary and relative). 
     198      !!       conserves the horizontal kinetic energy. 
     199      !!         The general trend of momentum is increased due to the vorticity  
     200      !!       term which is given by: 
     201      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ] 
     202      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ] 
     203      !!       where rvor is the relative vorticity at f-point 
     204      !! 
     205      !! ** Action : - Update (ua,va) with the now vorticity term trend 
     206      !!---------------------------------------------------------------------- 
     207      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index 
     208      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric 
     209      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities 
     210      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend 
     211      ! 
     212      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     213      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
     214      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zwt   ! 2D workspace 
     215      !!---------------------------------------------------------------------- 
     216      ! 
     217      IF( kt == nit000 ) THEN 
     218         IF(lwp) WRITE(numout,*) 
     219         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme' 
     220         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     221      ENDIF 
     222      ! 
     223      !                                                ! =============== 
     224      DO jk = 1, jpkm1                                 ! Horizontal slab 
     225         !                                             ! =============== 
     226         ! 
     227         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
     228         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     229            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk) 
     230         CASE ( np_RVO )                           !* relative vorticity 
     231            DO jj = 1, jpjm1 
     232               DO ji = 1, jpim1 
     233                  zwz(ji,jj) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
     234                     &          - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     235               END DO 
     236            END DO 
     237            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
     238               DO jj = 1, jpjm1 
     239                  DO ji = 1, jpim1 
     240                     zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     241                  END DO 
     242               END DO 
     243            ENDIF 
     244            CALL lbc_lnk( zwz, 'F', 1. ) 
     245            DO jj = 2, jpj 
     246               DO ji = 2, jpi   ! vector opt. 
     247                  zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ) + zwz(ji,jj  )   & 
     248                     &                  + zwz(ji-1,jj-1) + zwz(ji,jj-1)   ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     249               END DO 
     250            END DO 
     251         CASE ( np_MET )                           !* metric term 
     252            DO jj = 2, jpj 
     253               DO ji = 2, jpi 
     254                  zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
     255                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk) 
     256               END DO 
     257            END DO 
     258         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     259            DO jj = 1, jpjm1 
     260               DO ji = 1, jpim1                          ! relative vorticity 
     261                  zwz(ji,jj) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
     262                     &           - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
     263               END DO 
     264            END DO 
     265            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
     266               DO jj = 1, jpjm1 
     267                  DO ji = 1, jpim1 
     268                     zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     269                  END DO 
     270               END DO 
     271            ENDIF 
     272            CALL lbc_lnk( zwz, 'F', 1. ) 
     273            DO jj = 2, jpj 
     274               DO ji = 2, jpi   ! vector opt. 
     275                  zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ) + zwz(ji,jj  )    & 
     276                     &                                 + zwz(ji-1,jj-1) + zwz(ji,jj-1) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     277               END DO 
     278            END DO 
     279         CASE ( np_CME )                           !* Coriolis + metric 
     280            DO jj = 2, jpj 
     281               DO ji = 2, jpi   ! vector opt. 
     282                  zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
     283                       &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
     284                       &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk) 
     285               END DO 
     286            END DO 
     287         CASE DEFAULT                                             ! error 
     288            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
     289         END SELECT 
     290         ! 
     291         !                                   !==  compute and add the vorticity term trend  =! 
     292         DO jj = 2, jpjm1 
     293            DO ji = 2, jpim1   ! vector opt. 
     294               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    & 
     295                  &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   & 
     296                  &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   ) 
     297                  ! 
     298               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    & 
     299                  &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &  
     300                  &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )  
     301            END DO   
     302         END DO   
     303         !                                             ! =============== 
     304      END DO                                           !   End of slab 
     305      !                                                ! =============== 
     306   END SUBROUTINE vor_enT 
     307 
     308 
    166309   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 
    167310      !!---------------------------------------------------------------------- 
     
    217360            DO jj = 1, jpjm1 
    218361               DO ji = 1, fs_jpim1   ! vector opt. 
    219                   zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    220                        &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    221                        &     * 0.5 * r1_e1e2f(ji,jj) 
     362                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     363                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    222364               END DO 
    223365            END DO 
     
    225367            DO jj = 1, jpjm1 
    226368               DO ji = 1, fs_jpim1   ! vector opt. 
    227                   zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    228                      &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   & 
    229                      &                   * r1_e1e2f(ji,jj) 
     369                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      & 
     370                     &                        - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    230371               END DO 
    231372            END DO 
     
    233374            DO jj = 1, jpjm1 
    234375               DO ji = 1, fs_jpim1   ! vector opt. 
    235                   zwz(ji,jj) = ff_f(ji,jj)                                                                        & 
    236                        &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    237                        &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    238                        &     * 0.5 * r1_e1e2f(ji,jj) 
     376                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     377                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    239378               END DO 
    240379            END DO 
     
    328467            DO jj = 1, jpjm1 
    329468               DO ji = 1, fs_jpim1   ! vector opt. 
    330                   zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    331                        &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    332                        &     * 0.5 * r1_e1e2f(ji,jj) 
     469                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     470                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    333471               END DO 
    334472            END DO 
     
    336474            DO jj = 1, jpjm1 
    337475               DO ji = 1, fs_jpim1   ! vector opt. 
    338                   zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    339                      &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
    340                      &                   * r1_e1e2f(ji,jj) 
     476                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  & 
     477                     &                        - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    341478               END DO 
    342479            END DO 
     
    344481            DO jj = 1, jpjm1 
    345482               DO ji = 1, fs_jpim1   ! vector opt. 
    346                   zwz(ji,jj) = ff_f(ji,jj)                                                                       & 
    347                        &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    348                        &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    349                        &     * 0.5 * r1_e1e2f(ji,jj) 
     483                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     484                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    350485               END DO 
    351486            END DO 
     
    412547      INTEGER  ::   ierr         ! local integer 
    413548      REAL(wp) ::   zua, zva     ! local scalars 
    414       REAL(wp) ::   zmsk, ze3    ! local scalars 
    415       REAL(wp), DIMENSION(jpi,jpj)   :: zwx , zwy , zwz , z1_e3f 
    416       REAL(wp), DIMENSION(jpi,jpj)   :: ztnw, ztne, ztsw, ztse 
     549      REAL(wp) ::   zmsk, ze3f   ! local scalars 
     550      REAL(wp), DIMENSION(jpi,jpj) ::  zwx , zwy , zwz , z1_e3f 
     551      REAL(wp), DIMENSION(jpi,jpj) ::  ztnw, ztne, ztsw, ztse 
    417552      !!---------------------------------------------------------------------- 
    418553      ! 
     
    431566            DO jj = 1, jpjm1 
    432567               DO ji = 1, fs_jpim1   ! vector opt. 
    433                   ze3  = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     568                  ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    434569                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  ) 
    435                   IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3 
    436                   ELSE                      ;   z1_e3f(ji,jj) = 0._wp 
     570                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f 
     571                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp 
    437572                  ENDIF 
    438573               END DO 
     
    441576            DO jj = 1, jpjm1 
    442577               DO ji = 1, fs_jpim1   ! vector opt. 
    443                   ze3  = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     578                  ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    444579                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  ) 
    445580                  zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
    446581                     &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  ) 
    447                   IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3 
    448                   ELSE                      ;   z1_e3f(ji,jj) = 0._wp 
     582                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f 
     583                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp 
    449584                  ENDIF 
    450585               END DO 
     
    462597            DO jj = 1, jpjm1 
    463598               DO ji = 1, fs_jpim1   ! vector opt. 
    464                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    465                      &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
    466                      &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
     599                  zwz(ji,jj) = ( e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  & 
     600                     &         - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 
    467601               END DO 
    468602            END DO 
     
    470604            DO jj = 1, jpjm1 
    471605               DO ji = 1, fs_jpim1   ! vector opt. 
    472                   zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    473                        &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    474                        &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
     606                  zwz(ji,jj) = (   ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     607                     &           - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    475608               END DO 
    476609            END DO 
     
    478611            DO jj = 1, jpjm1 
    479612               DO ji = 1, fs_jpim1   ! vector opt. 
    480                   zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    481                      &                           - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
    482                      &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj) 
     613                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      & 
     614                     &                           - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )  & 
     615                     &                        * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    483616               END DO 
    484617            END DO 
     
    486619            DO jj = 1, jpjm1 
    487620               DO ji = 1, fs_jpim1   ! vector opt. 
    488                   zwz(ji,jj) = (  ff_f(ji,jj)                                                                        & 
    489                        &        + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    490                        &            - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    491                        &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
     621                  zwz(ji,jj) = (   ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     622                     &                         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    492623               END DO 
    493624            END DO 
     
    543674 
    544675 
     676 
     677   SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva ) 
     678      !!---------------------------------------------------------------------- 
     679      !!                ***  ROUTINE vor_eeT  *** 
     680      !! 
     681      !! ** Purpose :   Compute the now total vorticity trend and add it to  
     682      !!      the general trend of the momentum equation. 
     683      !! 
     684      !! ** Method  :   Trend evaluated using now fields (centered in time)  
     685      !!      and the Arakawa and Lamb (1980) vector form formulation using  
     686      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een). 
     687      !!      The change consists in  
     688      !!      Add this trend to the general momentum trend (ua,va). 
     689      !! 
     690      !! ** Action : - Update (ua,va) with the now vorticity term trend 
     691      !! 
     692      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
     693      !!---------------------------------------------------------------------- 
     694      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
     695      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
     696      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities 
     697      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend 
     698      ! 
     699      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     700      INTEGER  ::   ierr           ! local integer 
     701      REAL(wp) ::   zua, zva       ! local scalars 
     702      REAL(wp) ::   zmsk, z1_e3t   ! local scalars 
     703      REAL(wp), DIMENSION(jpi,jpj) ::   zwx , zwy , zwz 
     704      REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse 
     705      !!---------------------------------------------------------------------- 
     706      ! 
     707      IF( kt == nit000 ) THEN 
     708         IF(lwp) WRITE(numout,*) 
     709         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 
     710         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     711      ENDIF 
     712      ! 
     713      !                                                ! =============== 
     714      DO jk = 1, jpkm1                                 ! Horizontal slab 
     715         !                                             ! =============== 
     716         ! 
     717         ! 
     718         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     719         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     720            DO jj = 1, jpjm1 
     721               DO ji = 1, fs_jpim1   ! vector opt. 
     722                  zwz(ji,jj) = ff_f(ji,jj) 
     723               END DO 
     724            END DO 
     725         CASE ( np_RVO )                           !* relative vorticity 
     726            DO jj = 1, jpjm1 
     727               DO ji = 1, fs_jpim1   ! vector opt. 
     728                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     729                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
     730                     &       * r1_e1e2f(ji,jj) 
     731               END DO 
     732            END DO 
     733         CASE ( np_MET )                           !* metric term 
     734            DO jj = 1, jpjm1 
     735               DO ji = 1, fs_jpim1   ! vector opt. 
     736                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     737                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
     738               END DO 
     739            END DO 
     740         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     741            DO jj = 1, jpjm1 
     742               DO ji = 1, fs_jpim1   ! vector opt. 
     743                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     744                     &                           - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
     745                     &                      * r1_e1e2f(ji,jj)    ) 
     746               END DO 
     747            END DO 
     748         CASE ( np_CME )                           !* Coriolis + metric 
     749            DO jj = 1, jpjm1 
     750               DO ji = 1, fs_jpim1   ! vector opt. 
     751                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     752                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
     753               END DO 
     754            END DO 
     755         CASE DEFAULT                                             ! error 
     756            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
     757         END SELECT 
     758         ! 
     759         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
     760            DO jj = 1, jpjm1 
     761               DO ji = 1, fs_jpim1   ! vector opt. 
     762                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     763               END DO 
     764            END DO 
     765         ENDIF 
     766         ! 
     767         CALL lbc_lnk( zwz, 'F', 1. ) 
     768         ! 
     769         !                                   !==  horizontal fluxes  ==! 
     770         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
     771         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     772 
     773         !                                   !==  compute and add the vorticity term trend  =! 
     774         jj = 2 
     775         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0 
     776         DO ji = 2, jpi          ! split in 2 parts due to vector opt. 
     777               z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
     778               ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) * z1_e3t 
     779               ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) ) * z1_e3t 
     780               ztse(ji,jj) = ( zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t 
     781               ztsw(ji,jj) = ( zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) ) * z1_e3t 
     782         END DO 
     783         DO jj = 3, jpj 
     784            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3 
     785               z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
     786               ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) * z1_e3t 
     787               ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) ) * z1_e3t 
     788               ztse(ji,jj) = ( zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t 
     789               ztsw(ji,jj) = ( zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) ) * z1_e3t 
     790            END DO 
     791         END DO 
     792         DO jj = 2, jpjm1 
     793            DO ji = fs_2, fs_jpim1   ! vector opt. 
     794               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     795                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     796               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
     797                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     798               pua(ji,jj,jk) = pua(ji,jj,jk) + zua 
     799               pva(ji,jj,jk) = pva(ji,jj,jk) + zva 
     800            END DO   
     801         END DO   
     802         !                                             ! =============== 
     803      END DO                                           !   End of slab 
     804      !                                                ! =============== 
     805   END SUBROUTINE vor_eeT 
     806 
     807 
    545808   SUBROUTINE dyn_vor_init 
    546809      !!--------------------------------------------------------------------- 
     
    550813      !!              tracer advection schemes 
    551814      !!---------------------------------------------------------------------- 
    552       INTEGER ::   ioptio          ! local integer 
    553       INTEGER ::   ji, jj, jk      ! dummy loop indices 
    554       INTEGER ::   ios             ! Local integer output status for namelist read 
    555       !! 
    556       NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix,   & 
    557          &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_msk 
    558       !!---------------------------------------------------------------------- 
    559  
     815      INTEGER ::   ji, jj, jk    ! dummy loop indices 
     816      INTEGER ::   ioptio, ios   ! local integer 
     817      !! 
     818      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   & 
     819         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk 
     820      !!---------------------------------------------------------------------- 
     821      ! 
     822      IF(lwp) THEN 
     823         WRITE(numout,*) 
     824         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency' 
     825         WRITE(numout,*) '~~~~~~~~~~~~' 
     826      ENDIF 
     827      ! 
    560828      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options 
    561829      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901) 
     
    565833902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp ) 
    566834      IF(lwm) WRITE ( numond, namdyn_vor ) 
    567  
     835      ! 
    568836      IF(lwp) THEN                    ! Namelist print 
    569          WRITE(numout,*) 
    570          WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency' 
    571          WRITE(numout,*) '~~~~~~~~~~~~' 
    572837         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme' 
    573          WRITE(numout,*) '      energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene 
    574838         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens 
    575          WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix 
     839         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene 
     840         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT 
     841         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT 
    576842         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een 
    577843         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f 
     844         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix 
    578845         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
    579846      ENDIF 
     847 
     848      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available') 
    580849 
    581850!!gm  this should be removed when choosing a unique strategy for fmask at the coast 
     
    586855      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 
    587856         DO jk = 1, jpk 
    588             DO jj = 2, jpjm1 
    589                DO ji = 2, jpim1 
    590                   IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) & 
    591                       fmask(ji,jj,jk) = 1._wp 
     857            DO jj = 1, jpjm1 
     858               DO ji = 1, jpim1 
     859                  IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              & 
     860                     & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )  fmask(ji,jj,jk) = 1._wp 
    592861               END DO 
    593862            END DO 
    594863         END DO 
    595           ! 
    596           CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
    597           ! 
     864         ! 
     865         CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
     866         ! 
    598867      ENDIF 
    599868!!gm end 
    600869 
    601870      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme) 
    602       IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENE   ;   ENDIF 
    603       IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENS   ;   ENDIF 
    604       IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_MIX   ;   ENDIF 
    605       IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_EEN   ;   ENDIF 
     871      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF 
     872      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF 
     873      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF 
     874      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF 
     875      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF 
     876      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF 
    606877      ! 
    607878      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 
     
    622893         nrvm = np_MET        ! metric term 
    623894         ntot = np_CME        ! Coriolis + metric term 
     895         ! 
     896         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term: 
     897         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2 
     898            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) ) 
     899            DO jj = 2, jpjm1 
     900               DO ji = 2, jpim1 
     901                  di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp 
     902                  dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp 
     903               END DO 
     904            END DO 
     905            CALL lbc_lnk_multi( di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions 
     906            ! 
     907         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f) 
     908            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) ) 
     909            DO jj = 1, jpjm1 
     910               DO ji = 1, jpim1 
     911                  di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj) 
     912                  dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj) 
     913               END DO 
     914            END DO 
     915            CALL lbc_lnk_multi( di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. )   ! Lateral boundary conditions 
     916         END SELECT 
     917         ! 
    624918      END SELECT 
    625919       
     
    627921         WRITE(numout,*) 
    628922         SELECT CASE( nvor_scheme ) 
    629          CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme' 
    630          CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme' 
    631          CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme' 
    632          CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme' 
     923         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)' 
     924         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)' 
     925         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)' 
     926         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)' 
     927         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' 
     928         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)' 
    633929         END SELECT          
    634930      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.