Changeset 13154


Ignore:
Timestamp:
2020-06-24T15:31:32+02:00 (6 weeks ago)
Author:
gm
Message:

Alternative Directions on VOR and KE and bug fixed ln_vorlat

Location:
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90

    r13151 r13154  
    7474      CASE( np_VEC_c2  )      
    7575         CALL dyn_keg     ( kt, nn_dynkeg,      Kmm, puu, pvv, Krhs )    ! vector form : horizontal gradient of kinetic energy 
     76!!an SWE : w = 0 
    7677         CALL dyn_zad     ( kt,                 Kmm, puu, pvv, Krhs )    ! vector form : vertical advection 
    7778      CASE( np_FLX_c2  )  
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90

    r13151 r13154  
    2929 
    3030   PUBLIC   dyn_keg    ! routine called by step module 
     31   PUBLIC   dyn_kegAD   ! routine called by step module 
    3132    
    3233   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2     = 0   !: 2nd order centered scheme (standard scheme) 
     
    155156      ! 
    156157   END SUBROUTINE dyn_keg 
    157  
     158    
     159    
     160   SUBROUTINE dyn_kegAD( kt, kscheme, puu, pvv, pu_rhs, pv_rhs ) 
     161      !!---------------------------------------------------------------------- 
     162      !!                  ***  ROUTINE dyn_kegAD  *** 
     163      !! 
     164      !! ** Purpose :   Compute the now momentum trend due to the horizontal 
     165      !!      gradient of the horizontal kinetic energy and add it to the  
     166      !!      general momentum trend. 
     167      !! 
     168      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that  
     169      !!      conserve kinetic energy. Compute the now horizontal kinetic energy  
     170      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 
     171      !!              * kscheme = nkeg_HW : Hollingsworth correction following 
     172      !!      Arakawa (2001). The now horizontal kinetic energy is given by: 
     173      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((u(j+1)+u(j-1))/2)^2  ) 
     174      !!                    + mj-1(  2 * vn^2 + ((v(i+1)+v(i-1))/2)^2  ) ] 
     175      !!       
     176      !!      Take its horizontal gradient and add it to the general momentum 
     177      !!      trend. 
     178      !!         u(rhs) = u(rhs) - 1/e1u di[ zhke ] 
     179      !!         v(rhs) = v(rhs) - 1/e2v dj[ zhke ] 
     180      !! 
     181      !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend 
     182      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing 
     183      !! 
     184      !! ** References : Arakawa, A., International Geophysics 2001. 
     185      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 
     186      !!---------------------------------------------------------------------- 
     187      ! 
     188      INTEGER                                  , INTENT( in )  ::  kt               ! ocean time-step index 
     189      INTEGER                                  , INTENT( in )  ::  kscheme          ! =0/1/2   type of KEG scheme  
     190      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::  puu, pvv         ! ocean velocities at Kmm 
     191      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::  pu_rhs, pv_rhs   ! RHS  
     192      ! 
     193      INTEGER  ::   ji, jj, jk             ! dummy loop indices 
     194      REAL(wp) ::   zu, zv                   ! local scalars 
     195      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke 
     196      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv  
     197      !!---------------------------------------------------------------------- 
     198      ! 
     199      IF( ln_timing )   CALL timing_start('dyn_keg') 
     200      ! 
     201      IF( kt == nit000 ) THEN 
     202         IF(lwp) WRITE(numout,*) 
     203         IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme 
     204         IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
     205      ENDIF 
     206       
     207      zhke(:,:,jpk) = 0._wp 
     208 
     209      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
     210!!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2) 
     211      CASE ( nkeg_C2_wpo )                          !--  Standard scheme  --! 
     212         DO_3D_01_01( 1, jpkm1 ) 
     213            zu =  (   puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)   & 
     214               &    + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk)   ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) 
     215            zv =  (   pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)   & 
     216               &    + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk)   ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) 
     217            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 
     218         END_3D 
     219!!an45          
     220      ! 
     221      CASE ( nkeg_C2 )                          !--  Standard scheme  --! 
     222         DO_3D_01_01( 1, jpkm1 ) 
     223            zu =    puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)   & 
     224               &  + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk) 
     225            zv =    pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)   & 
     226               &  + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk) 
     227            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 
     228         END_3D 
     229!!an 00_00 ? 
     230      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --! 
     231         DO_3D_00_00( 1, jpkm1 ) 
     232            zu = 8._wp * ( puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)    & 
     233               &         + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk) )  & 
     234               &   +     ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) * ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) )   & 
     235               &   +     ( puu(ji  ,jj-1,jk) + puu(ji  ,jj+1,jk) ) * ( puu(ji  ,jj-1,jk) + puu(ji  ,jj+1,jk) ) 
     236               ! 
     237            zv = 8._wp * ( pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)    & 
     238               &         + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk) )  & 
     239               &  +      ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) * ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) )   & 
     240               &  +      ( pvv(ji-1,jj  ,jk) + pvv(ji+1,jj  ,jk) ) * ( pvv(ji-1,jj  ,jk) + pvv(ji+1,jj  ,jk) ) 
     241            zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 
     242         END_3D 
     243         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 
     244         ! 
     245      END SELECT  
     246      ! 
     247         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***! 
     248            ! 
     249            DO_3D_00_00( 1, jpkm1 ) 
     250               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 
     251               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 
     252            END_3D 
     253            ! 
     254         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***! 
     255            ! 
     256            DO_3D_00_00( 1, jpkm1 ) 
     257               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 
     258            END_3D 
     259            ! 
     260         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***! 
     261            ! 
     262            DO_3D_00_00( 1, jpkm1 ) 
     263               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 
     264            END_3D 
     265            ! 
     266         ENDIF 
     267      IF( ln_timing )   CALL timing_stop('dyn_kegAD') 
     268      ! 
     269   END SUBROUTINE dyn_kegAD 
    158270   !!====================================================================== 
    159271END MODULE dynkeg 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90

    r13151 r13154  
    2222   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation 
    2323   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity 
     24   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  alternate direction computation of vorticity tendancy 
     25   !!                 !                                 for ENS, ENE 
    2426   !!---------------------------------------------------------------------- 
    2527 
     
    4547   USE lib_mpp        ! MPP library 
    4648   USE timing         ! Timing 
     49!!anAD only 
     50   USE dynkeg, ONLY : dyn_kegAD 
     51   USE dynadv, ONLY : nn_dynkeg 
    4752 
    4853   IMPLICIT NONE 
     
    5358 
    5459   !                                   !!* Namelist namdyn_vor: vorticity term 
    55    LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS) 
    56    LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE) 
    57    LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT) 
    58    LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET) 
    59    LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN) 
    60    INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
    61    LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX) 
    62    LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
     60   LOGICAL, PUBLIC ::   ln_dynvor_ens      !: enstrophy conserving scheme          (ENS) 
     61   LOGICAL, PUBLIC ::   ln_dynvor_ens_adVO   = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad) 
     62   LOGICAL, PUBLIC ::   ln_dynvor_ens_adKE   = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad) 
     63   LOGICAL, PUBLIC ::   ln_dynvor_ens_adKEVO = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad) 
     64   LOGICAL, PUBLIC ::   ln_dynvor_ene      !: f-point energy conserving scheme     (ENE) 
     65   LOGICAL, PUBLIC ::   ln_dynvor_ene_adVO   = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad) 
     66   LOGICAL, PUBLIC ::   ln_dynvor_ene_adKE   = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad) 
     67   LOGICAL, PUBLIC ::   ln_dynvor_ene_adKEVO = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad) 
     68   LOGICAL, PUBLIC ::   ln_dynvor_enT      !: t-point energy conserving scheme     (ENT) 
     69   LOGICAL, PUBLIC ::   ln_dynvor_eeT      !: t-point energy conserving scheme     (EET) 
     70   LOGICAL, PUBLIC ::   ln_dynvor_een      !: energy & enstrophy conserving scheme (EEN) 
     71   INTEGER, PUBLIC ::      nn_een_e3f         !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
     72   LOGICAL, PUBLIC ::   ln_dynvor_mix      !: mixed scheme                         (MIX) 
     73   LOGICAL, PUBLIC ::   ln_dynvor_msk      !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
    6374 
    6475   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme 
     
    7081   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme 
    7182   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme 
    72  
    73    INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
     83!!an 
     84   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adKE   = 11   ! ENS scheme - AD scheme (KE only) 
     85   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adVO   = 12   ! ENS scheme - AD scheme (VOR only) 
     86   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adKEVO = 13   ! ENS scheme - AD scheme (KE+VOR) 
     87   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adKE   = 21   ! ENE scheme - AD scheme (KE only) 
     88   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adVO   = 22   ! ENE scheme - AD scheme (VOR only) 
     89   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adKEVO = 23   ! ENE scheme - AD scheme (KE+VOR) 
     90!!an       
     91    
     92!!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO 
     93   INTEGER, PUBLIC ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
    7494   !                               ! associated indices: 
    7595   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
     
    99119CONTAINS 
    100120 
    101    SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs ) 
     121   SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs) 
    102122      !!---------------------------------------------------------------------- 
    103123      !! 
     
    109129      !!               for futher diagnostics (l_trddyn=T) 
    110130      !!---------------------------------------------------------------------- 
    111       INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index 
    112       INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices 
    113       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation 
     131      INTEGER ::   ji, jj, jk   ! dummy loop indice 
     132      INTEGER                             , INTENT( in  ) ::   kt               ! ocean time-step index 
     133      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs, Kbb   ! ocean time level indices 
     134      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv         ! ocean velocity field and RHS of momentum equation 
    114135      ! 
    115136      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     137      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zuu, zvv 
    116138      !!---------------------------------------------------------------------- 
    117139      ! 
     
    167189            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
    168190         CASE( np_ENE )                        !* energy conserving scheme 
     191                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   
     192                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     193            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     194             
     195         CASE( np_ENE_adVO )                     !* energy conserving scheme with alternating direction 
     196            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     197                               
     198                              !==  Alternative direction - VOR only  ==! 
     199 
     200                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     201 
     202                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     203                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     204                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     205                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     206                             CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     207                             DEALLOCATE( zuu ) 
     208            ELSE 
     209                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     210 
     211                             ALLOCATE( zvv(jpi,jpj,jpk) )  
     212                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     213                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     214                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     215                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     216                             DEALLOCATE( zvv ) 
     217            ENDIF 
     218            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     219         CASE( np_ENE_adKE )                     !* energy conserving scheme with alternating direction 
     220            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     221                                  
     222                                 !==  Alternative direction - KEG only  ==! 
     223 
    169224                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     225 
     226                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     227                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     228                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     229                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     230                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     231                             DEALLOCATE( zuu ) 
     232            ELSE 
     233 
     234                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend                       
     235 
     236                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
     237                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     238                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     239                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     240                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm)  , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     241                             DEALLOCATE( zvv ) 
     242            ENDIF 
    170243            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     244             
     245         CASE( np_ENE_adKEVO )                     !* energy conserving scheme with alternating direction 
     246         ! equivalent ADE pas sur coriolis 
     247         ! appel vor_ene(ncor) 
     248            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     249                               
     250                              !==  Alternative direction - KE + VOR  ==! 
     251 
     252                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     253                             !puis call que sur nrvm  
     254                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     255                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   !  
     256                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     257                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     258                             !puis call que sur nrvm  
     259                             CALL vor_ene(   kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     260                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )    
     261                             DEALLOCATE( zuu ) 
     262            ELSE 
     263 
     264                             ALLOCATE( zvv(jpi,jpj,jpk) )  
     265                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     266                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )    
     267                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     268                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     269                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     270                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )    
     271                             DEALLOCATE( zvv ) 
     272            ENDIF    
    171273         CASE( np_ENS )                        !* enstrophy conserving scheme 
     274                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   
     275                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend 
     276            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend 
     277         CASE( np_ENS_adVO )                     !* enstrophy conserving scheme with alternating direction 
     278            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     279    
     280                                 !==  Alternative direction - VOR only  ==! 
     281 
     282                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   
     283 
     284                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     285                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     286                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     287                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     288                             CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     289                             DEALLOCATE( zuu ) 
     290            ELSE 
     291                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 
     292 
     293                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
     294                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     295                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     296                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     297                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     298                             DEALLOCATE( zvv ) 
     299            ENDIF 
     300            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend 
     301         CASE( np_ENS_adKE )                     !* enstrophy conserving scheme with alternating direction 
     302            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     303             
     304                                !==  Alternative direction - KEG only  ==! 
     305 
    172306                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend 
     307 
     308                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     309                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     310                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     311                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     312                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     313                             DEALLOCATE( zuu ) 
     314            ELSE 
     315 
     316                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend                                                           
     317 
     318                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
     319                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     320                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     321                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     322                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:)  , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     323                             DEALLOCATE( zvv ) 
     324            ENDIF 
     325         CASE( np_ENS_adKEVO )                     !* enstrophy conserving scheme with alternating direction 
     326            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     327                               
     328                              !==  Alternative direction - KE + VOR  ==! 
     329 
     330                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     331                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     332                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   !  
     333                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     334                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     335                             CALL vor_ens(   kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     336                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )    
     337                             DEALLOCATE( zuu ) 
     338            ELSE 
     339 
     340                             ALLOCATE( zvv(jpi,jpj,jpk) )  
     341                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     342                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )    
     343                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     344                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     345                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     346                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )    
     347                             DEALLOCATE( zvv ) 
     348            ENDIF    
    173349            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend 
    174350         CASE( np_MIX )                        !* mixed ene-ens scheme 
     
    319495      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    320496      !!---------------------------------------------------------------------- 
    321       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
    322       INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    323       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    324       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities 
    325       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
     497      INTEGER                                  , INTENT(in   ) ::   kt          ! ocean time-step index 
     498      INTEGER                                  , INTENT(in   ) ::   Kmm              ! ocean time level index 
     499      INTEGER                                  , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
     500      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::   pu, pv    ! now velocities 
     501      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
    326502      ! 
    327503      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    328504      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
    329505      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace 
     506      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t        ! unpenalised scale factors 
    330507      !!---------------------------------------------------------------------- 
    331508      ! 
     
    377554         END SELECT 
    378555         ! 
    379          !                                   !==  horizontal fluxes and potential vorticity ==! 
    380          zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    381          zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    382          zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    383          ! 
    384          !                                   !==  compute and add the vorticity term trend  =! 
    385          DO_2D_00_00 
    386             zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
    387             zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
    388             zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    389             zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    390             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    391             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
    392          END_2D 
     556         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***! 
     557            ! 
     558            !                                   !==  horizontal fluxes and potential vorticity ==! 
     559            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     560            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     561            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     562            ! 
     563            !                                   !==  compute and add the vorticity term trend  =! 
     564            DO_2D_00_00 
     565               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
     566               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
     567               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
     568               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
     569               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     570               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     571            END_2D 
     572            !             
     573            ! 
     574         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***! 
     575            ! 
     576            ! 
     577            !                                   !==  horizontal fluxes and potential vorticity ==! 
     578!!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée  
     579!      on transportera (en AD) de la bonne manière 
     580            !zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     581            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * vv(:,:,jk,Kmm) 
     582            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     583            ! 
     584            !                                   !==  compute and add the vorticity term trend  =! 
     585            DO_2D_00_00 
     586               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
     587               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
     588               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     589            END_2D 
     590            ! 
     591         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***! 
     592            ! 
     593            ! 
     594            !                                   !==  horizontal fluxes and potential vorticity ==! 
     595!!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée  
     596            !zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     597            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * uu(:,:,jk,Kmm) 
     598            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     599            ! 
     600            !                                   !==  compute and add the vorticity term trend  =! 
     601            DO_2D_00_00 
     602               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
     603               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
     604               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     605            END_2D 
     606            ! 
     607         ENDIF 
    393608         !                                             ! =============== 
    394609      END DO                                           !   End of slab 
     
    417632      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    418633      !!---------------------------------------------------------------------- 
    419       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
    420       INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    421       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    422       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities 
    423       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
     634      INTEGER                                  , INTENT(in   ) ::   kt          ! ocean time-step index 
     635      INTEGER                                  , INTENT(in   ) ::   Kmm              ! ocean time level index 
     636      INTEGER                                  , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
     637      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::   pu, pv    ! now velocities 
     638      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
    424639      ! 
    425640      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    475690         ! 
    476691         ! 
    477          !                                   !==  horizontal fluxes and potential vorticity ==! 
    478          zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    479          zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    480          zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    481          ! 
    482          !                                   !==  compute and add the vorticity term trend  =! 
    483          DO_2D_00_00 
    484             zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  & 
    485                &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  ) 
    486             zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  & 
    487                &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  ) 
    488             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
    489             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    490          END_2D 
     692!!an wut ? v et u  
     693         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***! 
     694            ! 
     695            !                                   !==  horizontal fluxes and potential vorticity ==! 
     696            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     697            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     698            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     699            ! 
     700            !                                   !==  compute and add the vorticity term trend  =! 
     701            DO_2D_00_00 
     702               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  & 
     703                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  ) 
     704               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  & 
     705                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  ) 
     706               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     707               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     708            END_2D 
     709            ! 
     710         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***! 
     711            ! 
     712            ! 
     713            !                                   !==  horizontal fluxes and potential vorticity ==! 
     714            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     715            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     716            ! 
     717            !                                   !==  compute and add the vorticity term trend  =! 
     718            DO_2D_00_00 
     719               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  & 
     720                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  ) 
     721               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     722            END_2D 
     723            ! 
     724         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***! 
     725            ! 
     726            ! 
     727            !                                   !==  horizontal fluxes and potential vorticity ==! 
     728            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     729            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     730            ! 
     731            !                                   !==  compute and add the vorticity term trend  =! 
     732            DO_2D_00_00 
     733               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  & 
     734                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  ) 
     735               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     736            END_2D 
     737            ! 
     738         ENDIF 
    491739         !                                             ! =============== 
    492740      END DO                                           !   End of slab 
     
    7721020      !! 
    7731021      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   & 
    774          &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk 
     1022         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk,   & 
     1023         &                 ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, &   ! Alternative direction parameters 
     1024         &                 ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO 
    7751025      !!---------------------------------------------------------------------- 
    7761026      ! 
     
    8091059         DO_3D_10_10( 1, jpk ) 
    8101060            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              & 
    811                & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp 
     1061               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp 
    8121062         END_3D 
    8131063         ! 
     
    8191069      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme) 
    8201070      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF 
     1071      IF( ln_dynvor_ens_adVO ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENS_adVO   ;   ENDIF 
     1072      IF( ln_dynvor_ens_adKE ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENS_adKE   ;   ENDIF 
     1073      IF( ln_dynvor_ens_adKEVO ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS_adKEVO   ;   ENDIF 
    8211074      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF 
     1075      IF( ln_dynvor_ene_adVO ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENE_adVO   ;   ENDIF 
     1076      IF( ln_dynvor_ene_adKE ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENE_adKE   ;   ENDIF 
     1077      IF( ln_dynvor_ene_adKEVO ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE_adKEVO   ;   ENDIF 
    8221078      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF 
    8231079      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF 
     
    8661122         WRITE(numout,*) 
    8671123         SELECT CASE( nvor_scheme ) 
    868          CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)' 
    869          CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)' 
    870          CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)' 
    871          CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)' 
    872          CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' 
    873          CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)' 
     1124          
     1125         CASE( np_ENS    )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)' 
     1126         CASE( np_ENS_adVO ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adVO) on vorticity only' 
     1127         CASE( np_ENS_adKE ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only' 
     1128         CASE( np_ENS_adKEVO ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity' 
     1129          
     1130         CASE( np_ENE    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)' 
     1131         CASE( np_ENE_adVO ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only' 
     1132         CASE( np_ENE_adKE ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only' 
     1133         CASE( np_ENE_adKEVO ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity' 
     1134          
     1135         CASE( np_ENT    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)' 
     1136         CASE( np_EET    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)' 
     1137         CASE( np_EEN    )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' 
     1138         CASE( np_MIX    )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)' 
    8741139         END SELECT          
    8751140      ENDIF 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90

    r13151 r13154  
    1717   USE phycst           ! physical constants 
    1818   USE usrdef_nam 
    19    ! 
     19   USE lib_mpp        ! MPP library 
    2020   USE iom              ! xIOs server  
    2121 
     
    6666      REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars 
    6767      REAL(wp) ::   zve3a, zve3n, zve3b    !   -      - 
    68       REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva 
     68      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva      
     69      REAL(wp), DIMENSION(jpi,jpj)    ::   zssh 
     70      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3u, ze3v 
    6971      !! --------------------------------------------------------------------- 
    7072#if defined key_agrif 
     
    141143               &            CALL Agrif_Sponge_dyn        ! momentum sponge 
    142144#endif 
    143                             CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
    144   
    145                             CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    146   
    147                             CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
    148145 
    149146!!an - calcul du gradient de pression horizontal (explicit) 
     147      zssh(:,:) = ssh(:,:,Nnn) 
     148# if defined key_bvp 
     149      !   gradient and divergence not penalised 
     150      zssh(:,:) = zssh(:,:) * r1_rpo(:,:,1)   
     151#endif  
    150152      DO_3D_00_00( 1, jpkm1 ) 
    151          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 
    152          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 
     153         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( zssh(ji+1,jj) - zssh(ji,jj) ) * r1_e1u(ji,jj) 
     154         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( zssh(ji,jj+1) - zssh(ji,jj) ) * r1_e2v(ji,jj) 
    153155      END_3D 
    154156      ! 
    155       ! add wind stress forcing and layer linear friction to the RHS  
     157       
     158!      IF( kstp == nit000 .AND. lwp ) THEN 
     159!         WRITE(numout,*) 
     160!         WRITE(numout,*) 'step.F90 : classic script used' 
     161!         WRITE(numout,*) '~~~~~~~' 
     162!         IF(       ln_dynvor_ens_adVO .OR. ln_dynvor_ens_adKE .OR. ln_dynvor_ens_adKEVO   & 
     163!         &    .OR. ln_dynvor_ene_adVO .OR. ln_dynvor_ene_adKE .OR. ln_dynvor_ene_adKEVO   ) THEN 
     164!            CALL ctl_stop('STOP','step : alternative direction asked but classis step'  ) 
     165!         ENDIF 
     166!      ENDIF 
     167!!an      
     168!                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)  ==> RHS 
     169!  
     170!                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! vorticity             ==> RHS 
     171!  
     172!!an     In dynvor, dynkegAD is called even if not AD, so we keep the same step.F90 
     173   
     174                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs)  ! vorticity            ==> RHS 
     175   
     176                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
     177 
     178      ! add wind stress forcing and layer linear friction to the RHS 
     179      ze3u(:,:,:) = e3u(:,:,:,Nnn) 
     180      ze3v(:,:,:) = e3v(:,:,:,Nnn) 
     181# if defined key_bvp 
     182      !ze3u(:,:,:) = ze3u(:,:,:) * r1_rpo(:,:,:)   
     183      !ze3v(:,:,:) = ze3v(:,:,:) * r1_rpo(:,:,:)   
     184#endif   
    156185      z1_2rho0 = 0.5_wp * r1_rho0 
    157186      DO_3D_00_00(1,jpkm1) 
    158          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   & 
     187         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / ze3u(ji,jj,jk)   & 
    159188            &                                  - rn_rfr * uu(ji,jj,jk,Nbb) 
    160          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   & 
    161             &                                  - rn_rfr * vv(ji,jj,jk,Nbb) 
     189         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ze3v(ji,jj,jk)   & 
     190            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)   
    162191      END_3D    
    163 !!an          
     192      ! 
     193# if defined key_bvp 
     194      !  Add frictionnal term   - sigma * u 
     195      ! can be done via sigmaT (mean_ij) 
     196      DO_3D_00_00(1,jpkm1) 
     197         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - bmpu(ji,jj,jk) * uu(ji,jj,jk,Nbb) 
     198         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - bmpv(ji,jj,jk) * vv(ji,jj,jk,Nbb) 
     199      END_3D  
     200#endif  
     201      !          
    164202   
    165203      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.