Ignore:
Timestamp:
2021-02-25T18:07:15+01:00 (5 months ago)
Author:
techene
Message:

#2605 RK3 time-stepping on OCE only (run on GYRE_PISCES without key_top)

Location:
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE
Files:
3 added
11 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domain.F90

    r14255 r14547  
    9191      INTEGER ::   ji, jj, jk, jt   ! dummy loop indices 
    9292      INTEGER ::   iconf = 0    ! local integers 
    93       REAL(wp)::   zrdt 
    9493      CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))" 
    9594      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level 
     
    372371      ! 
    373372      ! set current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3 
     373#if defined key_RK3 
     374      rDt   =         rn_Dt 
     375      r1_Dt = 1._wp / rDt 
     376      ! 
     377      IF(lwp) THEN 
     378         WRITE(numout,*) 
     379         WRITE(numout,*) '           ===>>>   Runge Kutta 3rd order (RK3) :   rDt = ', rDt 
     380         WRITE(numout,*) 
     381      ENDIF 
     382      ! 
     383#else 
    374384      rDt   = 2._wp * rn_Dt 
    375385      r1_Dt = 1._wp / rDt 
     386      ! 
     387      IF(lwp) THEN 
     388         WRITE(numout,*) 
     389         WRITE(numout,*) '           ===>>>   Modified Leap-Frog (MLF) :   rDt = ', rDt 
     390         WRITE(numout,*) 
     391      ENDIF 
     392      ! 
     393#endif 
    376394      ! 
    377395      IF( l_SAS .AND. .NOT.ln_linssh ) THEN 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/divhor.F90

    r13558 r14547  
    1313   !!             -   ! 2014-12  (G. Madec) suppression of cross land advection option 
    1414   !!             -   ! 2015-10  (G. Madec) add velocity and rnf flag in argument of div_hor 
     15   !!            4.5  ! 2015-10  (S. Techene, G. Madec)  hdiv replaced by e3divUh 
    1516   !!---------------------------------------------------------------------- 
    1617 
     
    3637   PRIVATE 
    3738 
    38    PUBLIC   div_hor    ! routine called by step.F90 and istate.F90 
     39   !                  !! * Interface 
     40   INTERFACE div_hor 
     41      MODULE PROCEDURE div_hor_RK3, div_hor_old 
     42   END INTERFACE 
     43 
     44   PUBLIC   div_hor    ! routine called by ssh_nxt.F90 and istate.F90 
    3945 
    4046   !! * Substitutions 
     
    4854CONTAINS 
    4955 
    50    SUBROUTINE div_hor( kt, Kbb, Kmm ) 
     56   SUBROUTINE div_hor_RK3( kt, Kbb, Kmm, puu, pvv, pe3divUh ) 
     57      !!---------------------------------------------------------------------- 
     58      !!                  ***  ROUTINE div_hor_RK3  *** 
     59      !!                     
     60      !! ** Purpose :   compute the horizontal divergence at now time-step 
     61      !! 
     62      !! ** Method  :   the now divergence is computed as : 
     63      !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
     64      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla)  
     65      !! 
     66      !! ** Action  : - thickness weighted horizontal divergence of in input velocity (puu,pvv) 
     67      !!---------------------------------------------------------------------- 
     68      INTEGER                         , INTENT(in   ) ::   kt, Kbb, Kmm   ! ocean time-step & time-level indices 
     69      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   puu, pvv       ! horizontal velocity 
     70      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pe3divUh       ! e3t*div[Uh] 
     71      ! 
     72      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     73      !!---------------------------------------------------------------------- 
     74      ! 
     75      IF( ln_timing )   CALL timing_start('div_hor_RK3') 
     76      ! 
     77      IF( kt == nit000 ) THEN 
     78         IF(lwp) WRITE(numout,*) 
     79         IF(lwp) WRITE(numout,*) 'div_hor_RK3 : thickness weighted horizontal divergence ' 
     80         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     81         hdiv    (:,:,:) = 0._wp    ! initialize hdiv & pe3divUh for the halos and jpk level at the first time step 
     82      ENDIF 
     83      !  
     84      pe3divUh(:,:,:) = 0._wp    !!gm to be applied to the halos only 
     85      ! 
     86      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                                   
     87         hdiv(ji,jj,jk) = (   e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * puu(ji  ,jj,jk)      & 
     88            &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk)      & 
     89            &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * pvv(ji,jj  ,jk)      & 
     90            &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * pvv(ji,jj-1,jk)  )   & 
     91            &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     92      END_3D 
     93      ! 
     94      IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )             !==  + runoffs divergence  ==! 
     95      ! 
     96#if defined key_asminc  
     97      IF( ln_sshinc .AND. ln_asmiau )   &                      !==  + SSH assimilation increment  ==! 
     98         &           CALL ssh_asm_div( kt, Kbb, Kmm, hdiv ) 
     99#endif 
     100      ! 
     101      IF( ln_isf )   CALL isf_hdiv( kt, Kmm, hdiv )            !==  + ice-shelf mass exchange ==! 
     102      ! 
     103!!gm Patch before suppression of hdiv from all modules that use it 
     104      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                            !==  e3t * Horizontal divergence  ==! 
     105         pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
     106      END_3D 
     107!!gm end 
     108      ! 
     109      CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp )   !   (no sign change) 
     110      ! 
     111      IF( ln_timing )   CALL timing_stop('div_hor_RK3') 
     112      ! 
     113   END SUBROUTINE div_hor_RK3 
     114 
     115 
     116   SUBROUTINE div_hor_old( kt, Kbb, Kmm ) 
    51117      !!---------------------------------------------------------------------- 
    52118      !!                  ***  ROUTINE div_hor  *** 
     
    64130      ! 
    65131      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    66       REAL(wp) ::   zraur, zdep   ! local scalars 
    67       REAL(wp), DIMENSION(jpi,jpj) :: ztmp 
    68132      !!---------------------------------------------------------------------- 
    69133      ! 
     
    98162      IF( ln_timing )   CALL timing_stop('div_hor') 
    99163      ! 
    100    END SUBROUTINE div_hor 
     164   END SUBROUTINE div_hor_old 
    101165    
    102166   !!====================================================================== 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg.F90

    r14225 r14547  
    5757CONTAINS 
    5858 
    59    SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV ) 
     59   SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) 
    6060      !!---------------------------------------------------------------------- 
    6161      !!                  ***  ROUTINE dyn_spg  *** 
     
    7979      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation 
    8080      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels 
    81       INTEGER , OPTIONAL                  , INTENT( in )  ::  k_only_ADV          ! only Advection in the RHS 
    8281      ! 
    8382      INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
    84       REAL(wp) ::   z2dt, zg_2, zintp, zgrho0r, zld   ! local scalars 
     83      REAL(wp) ::   zg_2, zintp, zgrho0r, zld    ! local scalars 
    8584      REAL(wp)             , DIMENSION(jpi,jpj) ::   zpgu, zpgv   ! 2D workspace 
    8685      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zpice 
     
    151150         IF( ln_wave .and. ln_bern_srfc ) THEN          !== Add J terms: depth-independent Bernoulli head 
    152151            DO_2D( 0, 0, 0, 0 ) 
    153                zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3] 
    154                zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj) 
     152               zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj)   !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3] 
     153               zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e2v(ji,jj) 
    155154            END_2D 
    156155         ENDIF 
     
    167166      SELECT CASE ( nspg )                   !== surface pressure gradient computed and add to the general trend ==! 
    168167      CASE ( np_EXP )   ;   CALL dyn_spg_exp( kt,      Kmm,       puu, pvv, Krhs )                    ! explicit 
    169       CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV ) ! time-splitting 
     168      CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting 
    170169      END SELECT 
    171170      ! 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg_ts.F90

    r14225 r14547  
    6868   PUBLIC dyn_spg_ts        ! called by dyn_spg  
    6969   PUBLIC dyn_spg_ts_init   !    -    - dyn_spg_init 
     70   PUBLIC dyn_drg_init      ! called by rk3ssh 
    7071 
    7172   !! Time filtered arrays at baroclinic time step: 
     
    8081   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftsw, ftse         ! (only used with een vorticity scheme) 
    8182 
     83   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshe_rhs           ! RHS of ssh Eq. 
     84   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ue_rhs, Ve_rhs    ! RHS of barotropic velocity Eq. 
     85   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   CdU_u, CdU_v       ! top/bottom stress at u- & v-points 
     86 
    8287   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! local ratios 
    8388   REAL(wp) ::   r1_8  = 0.125_wp         ! 
     
    117122 
    118123 
    119    SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV ) 
     124   SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) 
    120125      !!---------------------------------------------------------------------- 
    121126      !! 
     
    145150      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices 
    146151      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation 
    147       REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels 
    148       INTEGER , OPTIONAL                  , INTENT( in )  ::  k_only_ADV          ! only Advection in the RHS 
     152      REAL(wp), DIMENSION(jpi,jpj    ,jpt), INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels 
    149153      ! 
    150154      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
     
    152156      LOGICAL  ::   ll_init               ! =T : special startup of 2d equations 
    153157      INTEGER  ::   noffset               ! local integers  : time offset for bdy update 
    154       REAL(wp) ::   r1_Dt_b, z1_hu, z1_hv          ! local scalars 
    155       REAL(wp) ::   za0, za1, za2, za3              !   -      - 
     158      REAL(wp) ::   z1_hu, z1_hv              ! local scalars 
     159      REAL(wp) ::   za0, za1, za2, za3        !   -      - 
    156160      REAL(wp) ::   zztmp, zldg               !   -      - 
    157       REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      - 
    158       REAL(wp) ::   zun_save, zvn_save              !   -      - 
     161      REAL(wp) ::   zhu_bck, zhv_bck, zhdiv   !   -      - 
     162      REAL(wp) ::   zun_save, zvn_save        !   -      - 
    159163      REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc 
    160164      REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg 
     
    163167      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points 
    164168      REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes 
    165 !!st#if defined key_qco  
    166 !!st      REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 
    167 !!st#endif 
    168169      ! 
    169170      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.  
     
    175176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 
    176177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask 
     178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d          ! 2D workspace 
    177179      REAL(wp) ::   zt0substep !   Time of day at the beginning of the time substep 
    178180      !!---------------------------------------------------------------------- 
     
    185187!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 
    186188      !                                         ! inverse of baroclinic time step  
    187       r1_Dt_b = 1._wp / rDt  
    188189      ! 
    189190      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart  
     
    228229      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step) 
    229230      ! ----------------------------------------------------------------------------- 
     231#if defined key_RK3 
     232      !                    !========================================! 
     233      !                    !==  Phase 1 for RK3 time integration  ==! 
     234      !                    !========================================! 
     235      ! 
     236      !                          ! Currently, RK3 requires the forward mode and NO time filtering of barotropic variables 
     237      IF( kt == nit000 ) THEN 
     238         IF( .NOT.ln_bt_fw  .OR. ln_bt_av )   CALL ctl_stop( 'dyn_spg_ts: RK3 requires ln_bt_fw=T AND ln_bt_av=F') 
     239      ENDIF 
     240      ! 
     241      !                          ! set values computed in RK3_ssh 
     242      zssh_frc(:,:) = sshe_rhs(:,:) 
     243        zu_frc(:,:) =   Ue_rhs(:,:) 
     244        zv_frc(:,:) =   Ve_rhs(:,:) 
     245      zCdU_u  (:,:) = CdU_u   (:,:) 
     246      zCdU_v  (:,:) = CdU_v   (:,:) 
     247 
     248!!gm ==>>> !!ts    ISSUe her on en discute    changer les cas ENS ENE  et triad ? 
     249      IF( kt == nit000 .AND. .NOT. ln_linssh )   CALL dyn_cor_2D_init( Kmm )   ! Set zwz, the barotropic Coriolis force coefficient 
     250      !                      ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes 
     251      ! 
     252       
     253#else 
     254      !                    !========================================! 
     255      !                    !==  Phase 1 for MLF time integration  ==! 
     256      !                    !========================================! 
     257      ! 
    230258      !       
    231259      ! 
    232260      !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends) 
    233261      !                                   !  ---------------------------  ! 
    234 #if defined key_qco  
     262# if defined key_qco  
    235263      zu_frc(:,:) = SUM( e3u_0(:,:,:  ) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:) 
    236264      zv_frc(:,:) = SUM( e3v_0(:,:,:  ) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:) 
    237 #else 
     265# else 
    238266      zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu(:,:,Kmm) 
    239267      zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv(:,:,Kmm) 
    240 #endif  
     268# endif  
    241269      ! 
    242270      ! 
     
    256284      !                      ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes 
    257285      ! 
    258       IF( .NOT. PRESENT(k_only_ADV) ) THEN   !* remove the 2D Coriolis trend   
    259          zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes  
    260          zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column 
    261          ! 
    262          CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in 
    263             &                                                                          zu_trd, zv_trd   )   ! ==>> out 
    264          ! 
    265          DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term (and possibly spg) from barotropic trend 
    266             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 
    267             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 
    268          END_2D 
    269       ENDIF 
     286      zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes  
     287      zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column 
     288      ! 
     289      CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in 
     290         &                                                                          zu_trd, zv_trd   )   ! ==>> out 
     291      ! 
     292      DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term (and possibly spg) from barotropic trend 
     293         zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 
     294         zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 
     295      END_2D 
    270296      ! 
    271297      !                                   !=  Add bottom stress contribution from baroclinic velocities  =! 
    272298      !                                   !  -----------------------------------------------------------  ! 
    273       IF( PRESENT(k_only_ADV) ) THEN         !* only Advection in the RHS : provide the barotropic bottom drag coefficients 
    274          DO_2D( 0, 0, 0, 0 ) 
    275             zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 
    276             zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 
    277          END_2D 
    278       ELSE                !* remove baroclinic drag AND provide the barotropic drag coefficients 
    279          CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) 
    280       ENDIF 
     299      CALL dyn_drg_init( Kbb, Kmm, puu   , pvv   , puu_b , pvv_b  ,   &     !  <<= IN 
     300         &                         zu_frc, zv_frc, zCdU_u, zCdU_v )         !  =>> OUT 
    281301      ! 
    282302      !                                   !=  Add atmospheric pressure forcing  =! 
     
    348368      END IF 
    349369      ! 
    350 #if defined key_asminc 
     370# if defined key_asminc 
    351371      !                                   !=  Add the IAU weighted SSH increment  =! 
    352372      !                                   !  ------------------------------------  ! 
     
    354374         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 
    355375      ENDIF 
     376# endif 
     377 
     378      !                    !==  END of  Phase 1 for MLF time integration  ==! 
    356379#endif 
     380 
     381 
    357382      !                                   != Fill boundary data arrays for AGRIF 
    358383      !                                   ! ------------------------------------ 
     
    701726         sshb_e (:,:) = sshn_e(:,:) 
    702727         sshn_e (:,:) = ssha_e(:,:) 
    703  
    704          !                                             !* Sum over whole bt loop 
     728         ! 
     729         !                                             !* Sum over whole bt loop (except in weight average) 
    705730         !                                             !  ---------------------- 
    706          za1 = wgtbtp1(jn)                                     
    707          IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities 
    708             puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:)  
    709             pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:)  
    710          ELSE                                       ! Sum transports 
    711             IF ( .NOT.ln_wd_dl ) THEN   
    712                puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) 
    713                pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) 
    714             ELSE  
    715                puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:) 
    716                pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:) 
    717             END IF  
    718          ENDIF 
    719          !                                          ! Sum sea level 
    720          pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:) 
    721  
     731         IF( ln_bt_av ) THEN 
     732            za1 = wgtbtp1(jn)                                     
     733            IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities 
     734               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:)  
     735               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:)  
     736            ELSE                                       ! Sum transports 
     737               IF ( .NOT.ln_wd_dl ) THEN   
     738                  puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) 
     739                  pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) 
     740               ELSE  
     741                  puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:) 
     742                  pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:) 
     743               ENDIF 
     744            ENDIF 
     745            !                                          ! Sum sea level 
     746            pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:) 
     747         ENDIF 
     748         ! 
    722749         !                                                 ! ==================== ! 
    723750      END DO                                               !        end loop      ! 
    724751      !                                                    ! ==================== ! 
     752 
     753       
    725754      ! ----------------------------------------------------------------------------- 
    726755      ! Phase 3. update the general trend with the barotropic trend 
    727756      ! ----------------------------------------------------------------------------- 
    728757      ! 
    729       ! Set advection velocity correction: 
    730       IF (ln_bt_fw) THEN 
     758      IF(.NOT.ln_bt_av ) THEN                          !* Update Kaa barotropic external mode  
     759         puu_b(:,:,Kaa) = ua_e  (:,:) 
     760         pvv_b(:,:,Kaa) = va_e  (:,:) 
     761         pssh (:,:,Kaa) = ssha_e(:,:) 
     762      ENDIF 
     763 
     764#if defined key_RK3 
     765      !                                                !*  RK3 case 
     766      ! 
     767      ! advective transport from N to N+1 (i.e. Kbb to Kaa)  
     768      ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation (NOT USED) 
     769      vb2_b(:,:) = vn_adv(:,:) 
     770      !  
     771      IF( iom_use("ubar") ) THEN    ! RK3 single first: hu[N+1/2] = 1/2 ( hu[N] + hu[N+1] )  
     772         ALLOCATE( z2d(jpi,jpj) ) 
     773         z2d(:,:) = 2._wp / ( hu_e(:,:) + hu(:,:,Kbb) + 1._wp - ssumask(:,:) )  
     774         CALL iom_put(  "ubar", un_adv(:,:)*z2d(:,:) )    ! barotropic i-current 
     775         z2d(:,:) = 2._wp / ( hv_e(:,:) + hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 
     776         CALL iom_put(  "vbar", vn_adv(:,:)*z2d(:,:) )    ! barotropic i-current 
     777         DEALLOCATE( z2d ) 
     778      ENDIF 
     779      ! 
     780      !                    !==  END Phase 3 for RK3 (forward mode) ==! 
     781 
     782#else 
     783      !                                                !*  MLF case 
     784      ! 
     785      ! Set advective velocity correction: 
     786      IF( ln_bt_fw ) THEN 
    731787         IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 
    732788            DO_2D( 1, 1, 1, 1 ) 
     
    748804            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation 
    749805            vb2_b(:,:) = vn_adv(:,:) 
    750          END IF 
    751       ENDIF 
    752  
    753  
     806         ENDIF 
     807      ENDIF 
    754808      ! 
    755809      ! Update barotropic trend: 
    756810      IF( ln_dynadv_vec .OR. ln_linssh ) THEN 
    757811         DO jk=1,jpkm1 
    758             puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b 
    759             pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b 
     812            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt 
     813            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt 
    760814         END DO 
    761815      ELSE 
    762          ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 
    763 #if defined key_qcoTest_FluxForm 
    764          !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
    765          DO_2D( 1, 0, 1, 0 ) 
    766             zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj  ,Kaa) ) * ssumask(ji,jj) 
    767             zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji  ,jj+1,Kaa) ) * ssvmask(ji,jj) 
    768          END_2D 
    769 #else 
    770          DO_2D( 1, 0, 1, 0 ) 
    771             zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)   & 
    772                &                                      + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj) 
    773             zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)   & 
    774                &                                      + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 
    775          END_2D 
    776 #endif    
    777          CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
    778          ! 
    779          DO jk=1,jpkm1 
    780             puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   & 
    781                &             * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 
    782             pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   & 
    783                &             * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 
    784          END DO 
    785          ! Save barotropic velocities not transport: 
    786          puu_b(:,:,Kaa) =  puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
    787          pvv_b(:,:,Kaa) =  pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     816         IF(.NOT.ln_bt_av ) THEN   ! (puu_b,pvv_b)_Kaa is a velocity (hu,hv)_Kaa = (hu_e,hv_e) 
     817            !  
     818            DO jk=1,jpkm1 
     819               puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   & 
     820                  &             * ( puu_b(:,:,Kaa)*hu_e(:,:) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt 
     821               pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   & 
     822                  &             * ( pvv_b(:,:,Kaa)*hv_e(:,:) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt 
     823            END DO 
     824            ! 
     825         ELSE                  ! at this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 
     826            ! 
     827# if defined key_qcoTest_FluxForm 
     828            !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
     829            DO_2D( 1, 0, 1, 0 ) 
     830               zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj  ,Kaa) ) * ssumask(ji,jj) 
     831               zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji  ,jj+1,Kaa) ) * ssvmask(ji,jj) 
     832            END_2D 
     833# else 
     834            DO_2D( 1, 0, 1, 0 ) 
     835               zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)   & 
     836                  &                                      + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj) 
     837               zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)   & 
     838                  &                                      + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 
     839            END_2D 
     840# endif    
     841            CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
     842            ! 
     843            DO jk=1,jpkm1 
     844               puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   & 
     845                  &             * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt 
     846               pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   & 
     847                  &             * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt 
     848            END DO 
     849            ! Save barotropic velocities not transport: 
     850            puu_b(:,:,Kaa) =  puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
     851            pvv_b(:,:,Kaa) =  pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     852         ENDIF 
    788853      ENDIF 
    789854 
     
    795860      END DO 
    796861 
    797       IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN  
     862      IF( ln_wd_dl .AND. ln_wd_dl_bc) THEN  
    798863         DO jk = 1, jpkm1 
    799864            puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) & 
    800                        & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk)  
     865               &            + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk)  
    801866            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &  
    802                        & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)  
     867               &            + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)  
    803868         END DO 
    804       END IF  
    805  
    806        
     869      ENDIF 
     870             
    807871      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current 
    808872      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current 
     873 
     874      !                    !==  END Phase 3 for MLF time integration  ==! 
     875#endif 
     876 
    809877      ! 
    810878#if defined key_agrif 
     
    834902   END SUBROUTINE dyn_spg_ts 
    835903 
    836  
    837    SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) 
     904    
     905   SUBROUTINE ts_wgt( ll_av, ll_fw, Kpit, zwgt1, zwgt2) 
    838906      !!--------------------------------------------------------------------- 
    839907      !!                   ***  ROUTINE ts_wgt  *** 
     
    841909      !! ** Purpose : Set time-splitting weights for temporal averaging (or not) 
    842910      !!---------------------------------------------------------------------- 
    843       LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true. 
    844       LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true. 
    845       INTEGER, INTENT(inout) :: jpit      ! cycle length     
    846       REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, & ! Primary weights 
    847                                                          zwgt2    ! Secondary weights 
    848        
    849       INTEGER ::  jic, jn, ji                      ! temporary integers 
    850       REAL(wp) :: za1, za2 
    851       !!---------------------------------------------------------------------- 
    852  
     911      LOGICAL, INTENT(in   ) ::   ll_av     ! temporal averaging=.true. 
     912      LOGICAL, INTENT(in   ) ::   ll_fw     ! forward time splitting =.true. 
     913      INTEGER, INTENT(inout) ::   Kpit      ! cycle length 
     914      !! 
     915      INTEGER ::  jic, jn, ji   ! local integers 
     916      REAL(wp) :: za1, za2      ! loca scalars 
     917      REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, zwgt2   ! Primary & Secondary weights 
     918      !!---------------------------------------------------------------------- 
     919      ! 
    853920      zwgt1(:) = 0._wp 
    854921      zwgt2(:) = 0._wp 
    855  
    856       ! Set time index when averaged value is requested 
    857       IF (ll_fw) THEN  
    858          jic = nn_e 
    859       ELSE 
    860          jic = 2 * nn_e 
    861       ENDIF 
    862  
    863       ! Set primary weights: 
    864       IF (ll_av) THEN 
    865            ! Define simple boxcar window for primary weights  
    866            ! (width = nn_e, centered around jic)      
    867          SELECT CASE ( nn_bt_flt ) 
    868               CASE( 0 )  ! No averaging 
    869                  zwgt1(jic) = 1._wp 
    870                  jpit = jic 
    871  
    872               CASE( 1 )  ! Boxcar, width = nn_e 
    873                  DO jn = 1, 3*nn_e 
    874                     za1 = ABS(float(jn-jic))/float(nn_e)  
    875                     IF (za1 < 0.5_wp) THEN 
    876                       zwgt1(jn) = 1._wp 
    877                       jpit = jn 
    878                     ENDIF 
    879                  ENDDO 
    880  
    881               CASE( 2 )  ! Boxcar, width = 2 * nn_e 
    882                  DO jn = 1, 3*nn_e 
    883                     za1 = ABS(float(jn-jic))/float(nn_e)  
    884                     IF (za1 < 1._wp) THEN 
    885                       zwgt1(jn) = 1._wp 
    886                       jpit = jn 
    887                     ENDIF 
    888                  ENDDO 
    889               CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 
     922      ! 
     923      !                          !==  Set time index when averaged value is requested  ==! 
     924      IF (ll_fw) THEN   ;   jic =     nn_e 
     925      ELSE              ;   jic = 2 * nn_e 
     926      ENDIF 
     927      ! 
     928      !                          !==  Set primary weights  ==! 
     929      ! 
     930      IF (ll_av) THEN               != Define simple boxcar window for primary weights  
     931         !                                       ! (width = nn_e, centered around jic)      
     932         SELECT CASE( nn_bt_flt ) 
     933         ! 
     934         CASE( 0 )                  ! No averaging 
     935            zwgt1(jic) = 1._wp 
     936            Kpit = jic 
     937            ! 
     938         CASE( 1 )                  ! Boxcar, width = nn_e 
     939            DO jn = 1, 3*nn_e 
     940               za1 = ABS( REAL( jn-jic, wp) ) / REAL( nn_e, wp )  
     941               IF( za1 < 0.5_wp ) THEN 
     942                  zwgt1(jn) = 1._wp 
     943                  Kpit = jn 
     944               ENDIF 
     945            END DO 
     946            ! 
     947         CASE( 2 )                  ! Boxcar, width = 2 * nn_e 
     948            DO jn = 1, 3*nn_e 
     949               za1 = ABS(REAL( jn-jic, wp) ) / REAL( nn_e, wp )  
     950               IF( za1 < 1._wp ) THEN 
     951                  zwgt1(jn) = 1._wp 
     952                  Kpit = jn 
     953               ENDIF 
     954            END DO 
     955            ! 
     956         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 
     957         ! 
    890958         END SELECT 
    891959 
    892       ELSE ! No time averaging 
     960      ELSE                          != No time averaging 
    893961         zwgt1(jic) = 1._wp 
    894          jpit = jic 
    895       ENDIF 
    896      
    897       ! Set secondary weights 
    898       DO jn = 1, jpit 
    899         DO ji = jn, jpit 
    900              zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 
    901         END DO 
     962         Kpit = jic 
     963      ENDIF 
     964      ! 
     965      !                          !==  Set secondary weights  ==! 
     966      DO jn = 1, Kpit 
     967         DO ji = jn, Kpit 
     968            zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 
     969         END DO 
    902970      END DO 
    903  
    904       ! Normalize weigths: 
    905       za1 = 1._wp / SUM(zwgt1(1:jpit)) 
    906       za2 = 1._wp / SUM(zwgt2(1:jpit)) 
    907       DO jn = 1, jpit 
    908         zwgt1(jn) = zwgt1(jn) * za1 
    909         zwgt2(jn) = zwgt2(jn) * za2 
     971      ! 
     972      !                          !==  Normalize weigths  ==! 
     973      za1 = 1._wp / SUM( zwgt1(1:Kpit) ) 
     974      za2 = 1._wp / SUM( zwgt2(1:Kpit) ) 
     975      DO jn = 1, Kpit 
     976         zwgt1(jn) = zwgt1(jn) * za1 
     977         zwgt2(jn) = zwgt2(jn) * za2 
    910978      END DO 
    911979      ! 
     
    10231091         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_e in namelist   nn_e = ', nn_e 
    10241092      ENDIF 
    1025  
     1093      ! 
    10261094      IF(ln_bt_av) THEN 
    10271095         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_e time steps is on ' 
     
    10291097         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables ' 
    10301098      ENDIF 
    1031       ! 
    10321099      ! 
    10331100      IF(ln_bt_fw) THEN 
     
    10881155      !! and update depths at T- points (ht) at each barotropic time step 
    10891156      !! 
    1090       !! Compute zwz = f / ( height of the water colomn ) 
     1157      !! Compute zwz = f/h              (potential planetary voricity) 
     1158      !! Compute ftne, ftnw, ftse, ftsw (triad of potential planetary voricity) 
    10911159      !!---------------------------------------------------------------------- 
    10921160      INTEGER,  INTENT(in)         ::  Kmm  ! Time index 
     
    11381206         ! 
    11391207      END SELECT 
    1140        
     1208      ! 
    11411209   END SUBROUTINE dyn_cor_2d_init 
    11421210 
     
    12871355      !! ** Purpose :  
    12881356      !!---------------------------------------------------------------------- 
     1357      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn 
     1358      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 
     1359      ! 
    12891360      INTEGER  ::   ji ,jj               ! dummy loop indices 
    12901361      LOGICAL  ::   ll_tmp1, ll_tmp2 
    1291       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn 
    1292       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 
    1293       !!---------------------------------------------------------------------- 
     1362      !!---------------------------------------------------------------------- 
     1363      ! 
    12941364      DO_2D( 0, 0, 0, 0 ) 
    12951365         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                & 
     
    13301400         ENDIF 
    13311401      END_2D 
    1332              
     1402      ! 
    13331403   END SUBROUTINE wad_spg 
    13341404      
     
    13741444      ! 
    13751445      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities 
    1376           
    13771446         DO_2D( 0, 0, 0, 0 ) 
    13781447            ikbu = mbku(ji,jj)        
     
    13811450            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 
    13821451         END_2D 
    1383       ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities 
    1384           
     1452      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities         
    13851453         DO_2D( 0, 0, 0, 0 ) 
    13861454            ikbu = mbku(ji,jj)        
     
    14001468         END_2D 
    14011469      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 
    1402           
    14031470         DO_2D( 0, 0, 0, 0 ) 
    14041471            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 
     
    14121479         ! 
    14131480         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity 
    1414              
    14151481            DO_2D( 0, 0, 0, 0 ) 
    14161482               iktu = miku(ji,jj) 
     
    14191485               zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 
    14201486            END_2D 
    1421          ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity 
    1422              
     1487         ELSE                               ! CENTRED integration: use BEFORE top baroclinic velocity 
    14231488            DO_2D( 0, 0, 0, 0 ) 
    14241489               iktu = miku(ji,jj) 
     
    14281493            END_2D 
    14291494         ENDIF 
    1430          ! 
    1431          !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 
    1432           
     1495         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)        
    14331496         DO_2D( 0, 0, 0, 0 ) 
    14341497            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 
     
    14391502      ! 
    14401503   END SUBROUTINE dyn_drg_init 
     1504 
    14411505 
    14421506   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in 
     
    14771541   END SUBROUTINE ts_bck_interp 
    14781542 
    1479  
    14801543   !!====================================================================== 
    14811544END MODULE dynspg_ts 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90

    r14233 r14547  
    5151   IMPLICIT NONE 
    5252   PRIVATE 
     53    
     54   INTERFACE dyn_vor 
     55      MODULE PROCEDURE dyn_vor_3D, dyn_vor_RK3 
     56   END INTERFACE 
    5357 
    5458   PUBLIC   dyn_vor        ! routine called by step.F90 
     
    7478   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme 
    7579 
    76    INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
    77    !                               ! associated indices: 
     80   !                                    !: choice of calculated vorticity 
     81   INTEGER, PUBLIC ::   ncor, nrvm, ntot   ! Coriolis, relative vorticity, total vorticity 
     82   !                                       ! associated indices: 
    7883   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
    7984   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity 
     
    103108   !!---------------------------------------------------------------------- 
    104109CONTAINS 
    105  
    106    SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs ) 
     110    
     111   SUBROUTINE dyn_vor_3D( kt, Kmm, puu, pvv, Krhs ) 
    107112      !!---------------------------------------------------------------------- 
    108113      !! 
     
    121126      !!---------------------------------------------------------------------- 
    122127      ! 
    123       IF( ln_timing )   CALL timing_start('dyn_vor') 
     128      IF( ln_timing )   CALL timing_start('dyn_vor_3D') 
    124129      ! 
    125130      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==! 
     
    209214         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    210215      ! 
    211       IF( ln_timing )   CALL timing_stop('dyn_vor') 
    212       ! 
    213    END SUBROUTINE dyn_vor 
     216      IF( ln_timing )   CALL timing_stop('dyn_vor_3D') 
     217      ! 
     218   END SUBROUTINE dyn_vor_3D 
     219 
     220 
     221   SUBROUTINE dyn_vor_RK3( kt, Kmm, puu, pvv, Krhs, knoco ) 
     222      !!---------------------------------------------------------------------- 
     223      !! 
     224      !! ** Purpose :   compute the lateral ocean tracer physics. 
     225      !! 
     226      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 
     227      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
     228      !!               and planetary vorticity trends) and send them to trd_dyn 
     229      !!               for futher diagnostics (l_trddyn=T) 
     230      !!---------------------------------------------------------------------- 
     231      INTEGER                             , INTENT(in   ) ::   kt          ! ocean time-step index 
     232      INTEGER                             , INTENT(in   ) ::   Kmm, Krhs   ! ocean time level indices 
     233      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation 
     234      INTEGER                             , INTENT(in   ) ::   knoco       ! specified vorticity trend (= np_MET or np_RVO) 
     235      !!---------------------------------------------------------------------- 
     236      ! 
     237      IF( ln_timing )   CALL timing_start('dyn_vor_RK3') 
     238      ! 
     239      !              !==  total vorticity trend added to the general trend  ==! 
     240      !!st WARNING 22/02 !!!!!!!! stoke drift or not stoke drift ? => bar to do later !!! 
     241      !! stoke drift a garder pas vortex force a priori !! 
     242      !! ATTENTION déja appelé avec Kbb !! 
     243       
     244         ! 
     245         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==! 
     246         CASE( np_ENT )                        !* energy conserving scheme  (T-pts) 
     247                             CALL vor_enT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     248            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN 
     249                             CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     250            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN 
     251                             CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force 
     252            ENDIF 
     253         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t) 
     254                             CALL vor_eeT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     255            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN 
     256                             CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     257            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN 
     258                             CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force 
     259            ENDIF 
     260         CASE( np_ENE )                        !* energy conserving scheme 
     261                             CALL vor_ene( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     262            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN 
     263                             CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     264            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN 
     265                             CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force 
     266            ENDIF 
     267         CASE( np_ENS )                        !* enstrophy conserving scheme 
     268                             CALL vor_ens( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend 
     269 
     270            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN 
     271                             CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend 
     272            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN 
     273                             CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend and vortex force 
     274            ENDIF 
     275         CASE( np_MIX )                        !* mixed ene-ens scheme 
     276                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens) 
     277            IF( ln_stcor )        CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )        ! add the Stokes-Coriolis trend 
     278            IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add vortex force 
     279         CASE( np_EEN )                        !* energy and enstrophy conserving scheme 
     280                             CALL vor_een( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     281            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN 
     282                             CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     283            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN 
     284                             CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force 
     285            ENDIF 
     286         END SELECT 
     287      ! 
     288      !                       ! print sum trends (used for debugging) 
     289      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               & 
     290         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     291      ! 
     292      IF( ln_timing )   CALL timing_stop('dyn_vor_RK3') 
     293      ! 
     294   END SUBROUTINE dyn_vor_RK3 
    214295 
    215296 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynzdf.F90

    r13497 r14547  
    7070      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation 
    7171      ! 
    72       INTEGER  ::   ji, jj, jk         ! dummy loop indices 
    73       INTEGER  ::   iku, ikv           ! local integers 
    74       REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars 
    75       REAL(wp) ::   zzws, ze3va        !   -      - 
    76       REAL(wp) ::   z1_e3ua, z1_e3va   !   -      - 
    77       REAL(wp) ::   zWu , zWv          !   -      - 
    78       REAL(wp) ::   zWui, zWvi         !   -      - 
    79       REAL(wp) ::   zWus, zWvs         !   -      - 
     72      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     73      INTEGER  ::   iku, ikv             ! local integers 
     74      REAL(wp) ::   zzwi, ze3ua, zDt_2   ! local scalars 
     75      REAL(wp) ::   zzws, ze3va          !   -      - 
     76      REAL(wp) ::   z1_e3ua, z1_e3va     !   -      - 
     77      REAL(wp) ::   zWu , zWv            !   -      - 
     78      REAL(wp) ::   zWui, zWvi           !   -      - 
     79      REAL(wp) ::   zWus, zWvs           !   -      - 
    8080      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace  
    8181      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      - 
     
    9393         ENDIF 
    9494      ENDIF 
     95      ! 
     96      zDt_2 = rDt * 0.5_wp 
     97      ! 
    9598      !                             !* explicit top/bottom drag case 
    9699      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa)) 
     
    138141            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    & 
    139142               &             + r_vvl   * e3v(ji,jj,ikv,Kaa) 
    140             puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 
    141             pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
     143            puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 
     144            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
    142145         END_2D 
    143146         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF) 
     
    149152               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    & 
    150153                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa) 
    151                puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 
    152                pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
     154               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 
     155               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
    153156            END_2D 
    154157         END IF 
     
    158161      ! 
    159162      !                    !* Matrix construction 
    160       zdt = rDt * 0.5 
    161163      IF( ln_zad_Aimp ) THEN   !! 
    162164         SELECT CASE( nldf_dyn ) 
     
    165167               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    & 
    166168                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
    167                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     169               zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    168170                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
    169                zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     171               zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    170172                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
    171173               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
    172174               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
    173                zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )  
    174                zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
    175                zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     175               zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )  
     176               zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp ) 
     177               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
    176178            END_3D 
    177179         CASE DEFAULT               ! iso-level lateral mixing 
     
    179181               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &    ! after scale factor at U-point 
    180182                  &             + r_vvl   * e3u(ji,jj,jk,Kaa) 
    181                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   & 
     183               zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   & 
    182184                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
    183                zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   & 
     185               zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   & 
    184186                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
    185187               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua 
    186188               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 
    187                zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
    188                zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 
    189                zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
     189               zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp ) 
     190               zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp ) 
     191               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 
    190192            END_3D 
    191193         END SELECT 
     
    194196            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    & 
    195197               &             + r_vvl   * e3u(ji,jj,1,Kaa) 
    196             zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   & 
     198            zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   & 
    197199               &         / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 
    198200            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua 
    199             zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 
    200             zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 
     201            zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWus, 0._wp ) 
     202            zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) ) 
    201203         END_2D 
    202204      ELSE 
     
    206208               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    & 
    207209                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
    208                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     210               zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    209211                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
    210                zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     212               zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    211213                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
    212214               zwi(ji,jj,jk) = zzwi 
     
    218220               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    & 
    219221                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
    220                zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    & 
     222               zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    & 
    221223                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
    222                zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    & 
     224               zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    & 
    223225                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
    224226               zwi(ji,jj,jk) = zzwi 
     
    245247            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    & 
    246248               &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point 
    247             zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
     249            zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
    248250         END_2D 
    249251         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit) 
     
    253255               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    & 
    254256                  &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point 
    255                zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
     257               zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
    256258            END_2D 
    257259         END IF 
     
    280282         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    & 
    281283            &             + r_vvl   * e3u(ji,jj,1,Kaa)  
    282          puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     284         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    283285            &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1)  
    284286      END_2D 
     
    297299      ! 
    298300      !                       !* Matrix construction 
    299       zdt = rDt * 0.5 
    300301      IF( ln_zad_Aimp ) THEN   !! 
    301302         SELECT CASE( nldf_dyn ) 
     
    304305               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    & 
    305306                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
    306                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     307               zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    307308                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
    308                zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     309               zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    309310                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
    310311               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
    311312               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
    312                zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 
    313                zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 
    314                zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     313               zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp ) 
     314               zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp ) 
     315               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
    315316            END_3D 
    316317         CASE DEFAULT               ! iso-level lateral mixing 
     
    318319               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    & 
    319320                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
    320                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    & 
     321               zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    & 
    321322                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
    322                zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    & 
     323               zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    & 
    323324                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
    324325               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va 
    325326               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 
    326                zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp ) 
    327                zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp ) 
    328                zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
     327               zwi(ji,jj,jk) = zzwi  + zDt_2 * MIN( zWvi, 0._wp ) 
     328               zws(ji,jj,jk) = zzws  - zDt_2 * MAX( zWvs, 0._wp ) 
     329               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 
    329330            END_3D 
    330331         END SELECT 
     
    333334            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    & 
    334335               &             + r_vvl   * e3v(ji,jj,1,Kaa) 
    335             zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    & 
     336            zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    & 
    336337               &         / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 
    337338            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va 
    338             zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 
    339             zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 
     339            zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp ) 
     340            zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) ) 
    340341         END_2D 
    341342      ELSE 
     
    345346               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    & 
    346347                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
    347                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     348               zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    348349                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
    349                zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     350               zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    350351                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
    351352               zwi(ji,jj,jk) = zzwi 
     
    357358               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    & 
    358359                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
    359                zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    & 
     360               zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    & 
    360361                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
    361                zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    & 
     362               zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    & 
    362363                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
    363364               zwi(ji,jj,jk) = zzwi 
     
    383384            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    & 
    384385               &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point 
    385             zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
     386            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
    386387         END_2D 
    387388         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN 
     
    390391               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    & 
    391392                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point 
    392                zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 
     393               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 
    393394            END_2D 
    394395         ENDIF 
     
    417418         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    & 
    418419            &             + r_vvl   * e3v(ji,jj,1,Kaa)  
    419          pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     420         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    420421            &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1)  
    421422      END_2D 
     
    432433      ! 
    433434      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    434          ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:) 
    435          ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:) 
     435         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:) 
     436         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) )*r1_Dt - ztrdv(:,:,:) 
    436437         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm ) 
    437438         DEALLOCATE( ztrdu, ztrdv )  
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/sshwzv.F90

    r14205 r14547  
    1717   !!   ssh_nxt       : after ssh 
    1818   !!   ssh_atf       : time filter the ssh arrays 
    19    !!   wzv           : compute now vertical velocity 
     19   !!   wzv           : generic interface of vertical velocity calculation 
     20   !!   wzv_MLF       : MLF: compute NOW vertical velocity 
     21   !!   wzv_RK3       : RK3: Compute a vertical velocity 
    2022   !!---------------------------------------------------------------------- 
    2123   USE oce            ! ocean dynamics and tracers variables 
     
    4446   IMPLICIT NONE 
    4547   PRIVATE 
     48   !                  !! * Interface 
     49   INTERFACE wzv 
     50      MODULE PROCEDURE wzv_MLF, wzv_RK3 
     51   END INTERFACE 
    4652 
    4753   PUBLIC   ssh_nxt        ! called by step.F90 
     
    7783      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index 
    7884      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height 
    79       !  
     85      ! 
    8086      INTEGER  ::   jk      ! dummy loop index 
    8187      REAL(wp) ::   zcoef   ! local scalar 
    8288      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace 
     89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace 
    8390      !!---------------------------------------------------------------------- 
    8491      ! 
     
    9299      ! 
    93100      zcoef = 0.5_wp * r1_rho0 
    94  
    95101      !                                           !------------------------------! 
    96102      !                                           !   After Sea Surface Height   ! 
    97103      !                                           !------------------------------! 
    98104      IF(ln_wd_il) THEN 
    99          CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv ) 
     105         CALL wad_lmt( pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv ) 
    100106      ENDIF 
    101107 
     
    133139   END SUBROUTINE ssh_nxt 
    134140 
    135     
    136    SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww ) 
    137       !!---------------------------------------------------------------------- 
    138       !!                ***  ROUTINE wzv  *** 
     141 
     142   SUBROUTINE wzv_MLF( kt, Kbb, Kmm, Kaa, pww ) 
     143      !!---------------------------------------------------------------------- 
     144      !!                ***  ROUTINE wzv_MLF  *** 
    139145      !!                    
    140146      !! ** Purpose :   compute the now vertical velocity 
     
    157163      !!---------------------------------------------------------------------- 
    158164      ! 
    159       IF( ln_timing )   CALL timing_start('wzv') 
     165      IF( ln_timing )   CALL timing_start('wzv_MLF') 
    160166      ! 
    161167      IF( kt == nit000 ) THEN 
    162168         IF(lwp) WRITE(numout,*) 
    163          IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 
    164          IF(lwp) WRITE(numout,*) '~~~~~ ' 
     169         IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity ' 
     170         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    165171         ! 
    166172         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
     
    193199         END DO 
    194200         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 
    195          DEALLOCATE( zhdiv )  
     201         DEALLOCATE( zhdiv ) 
    196202         !                                            !=================================! 
    197203      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==! 
     
    204210         !                                            !==========================================! 
    205211         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence 
    206             pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 & 
    207                &                            + r1_Dt * (  e3t(:,:,jk,Kaa)        & 
    208                &                                       - e3t(:,:,jk,Kbb)  )   ) * tmask(:,:,jk) 
     212            !                                               ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]] 
     213!!gm old 
     214!            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)   & 
     215!               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)      & 
     216!               &                                       - e3t(:,:,jk,Kbb)  )   ) * tmask(:,:,jk) 
     217!!gm slightly faster : 
     218            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                            & 
     219               &                            + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) )  ) * tmask(:,:,jk) 
     220!!gm end 
     221 
     222 
    209223         END DO 
    210224      ENDIF 
     
    256270#endif 
    257271      ! 
    258       IF( ln_timing )   CALL timing_stop('wzv') 
    259       ! 
    260    END SUBROUTINE wzv 
     272      IF( ln_timing )   CALL timing_stop('wzv_MLF') 
     273      ! 
     274   END SUBROUTINE wzv_MLF 
     275 
     276 
     277   SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww ) 
     278      !!---------------------------------------------------------------------- 
     279      !!                ***  ROUTINE wzv_RK3  *** 
     280      !!                    
     281      !! ** Purpose :   compute the now vertical velocity 
     282      !! 
     283      !! ** Method  : - Using the incompressibility hypothesis, the vertical  
     284      !!      velocity is computed by integrating the horizontal divergence   
     285      !!      from the bottom to the surface minus the scale factor evolution. 
     286      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     287      !! 
     288      !! ** action  :   pww      : now vertical velocity 
     289      !! 
     290      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     291      !!---------------------------------------------------------------------- 
     292      INTEGER                         , INTENT(in   ) ::   kt             ! time step 
     293      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level indices 
     294      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   puu, pvv       ! horizontal velocity at Kmm 
     295      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::             pww  ! vertical velocity at Kmm 
     296      ! 
     297      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     298      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv 
     299      REAL(wp)             , DIMENSION(jpi,jpj,jpk) ::   ze3div 
     300      !!---------------------------------------------------------------------- 
     301      ! 
     302      IF( ln_timing )   CALL timing_start('wzv_RK3') 
     303      ! 
     304      IF( kt == nit000 ) THEN 
     305         IF(lwp) WRITE(numout,*) 
     306         IF(lwp) WRITE(numout,*) 'wzv_RK3 : now vertical velocity ' 
     307         IF(lwp) WRITE(numout,*) '~~~~~ ' 
     308         ! 
     309         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
     310      ENDIF 
     311      ! 
     312      CALL div_hor( kt, Kbb, Kmm, puu, pvv, ze3div ) 
     313      !                                           !------------------------------! 
     314      !                                           !     Now Vertical Velocity    ! 
     315      !                                           !------------------------------! 
     316      ! 
     317      !                                               !===============================! 
     318      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      !==  z_tilde and layer cases  ==! 
     319         !                                            !===============================! 
     320         ALLOCATE( zhdiv(jpi,jpj,jpk) )  
     321         ! 
     322         DO jk = 1, jpkm1 
     323            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
     324            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
     325            DO_2D( 0, 0, 0, 0 ) 
     326               zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
     327            END_2D 
     328         END DO 
     329         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
     330         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells? 
     331         !                             ! Same question holds for hdiv. Perhaps just for security 
     332         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     333            ! computation of w 
     334            pww(:,:,jk) = pww(:,:,jk+1) - (   ze3div(:,:,jk) + zhdiv(:,:,jk)   & 
     335               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)       & 
     336               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk) 
     337         END DO 
     338         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 
     339         DEALLOCATE( zhdiv )  
     340         !                                            !=================================! 
     341      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==! 
     342         !                                            !=================================! 
     343         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence 
     344            pww(:,:,jk) = pww(:,:,jk+1) - ze3div(:,:,jk) 
     345         END DO 
     346         !                                            !==========================================! 
     347      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco') 
     348         !                                            !==========================================! 
     349         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence 
     350            !                                               ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]] 
     351            pww(:,:,jk) = pww(:,:,jk+1) - (  ze3div(:,:,jk)                            & 
     352               &                            + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) )  ) * tmask(:,:,jk) 
     353         END DO 
     354      ENDIF 
     355 
     356      IF( ln_bdy ) THEN 
     357         DO jk = 1, jpkm1 
     358            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:) 
     359         END DO 
     360      ENDIF 
     361      ! 
     362#if defined key_agrif 
     363      IF( .NOT. AGRIF_Root() ) THEN 
     364         ! 
     365         ! Mask vertical velocity at first/last columns/row  
     366         ! inside computational domain (cosmetic)  
     367         DO jk = 1, jpkm1 
     368            IF( lk_west ) THEN                             ! --- West --- ! 
     369               DO ji = mi0(2+nn_hls), mi1(2+nn_hls) 
     370                  DO jj = 1, jpj 
     371                     pww(ji,jj,jk) = 0._wp  
     372                  END DO 
     373               END DO 
     374            ENDIF 
     375            IF( lk_east ) THEN                             ! --- East --- ! 
     376               DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls) 
     377                  DO jj = 1, jpj 
     378                     pww(ji,jj,jk) = 0._wp 
     379                  END DO 
     380               END DO 
     381            ENDIF 
     382            IF( lk_south ) THEN                            ! --- South --- ! 
     383               DO jj = mj0(2+nn_hls), mj1(2+nn_hls) 
     384                  DO ji = 1, jpi 
     385                     pww(ji,jj,jk) = 0._wp 
     386                  END DO 
     387               END DO 
     388            ENDIF 
     389            IF( lk_north ) THEN                            ! --- North --- ! 
     390               DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls) 
     391                  DO ji = 1, jpi 
     392                     pww(ji,jj,jk) = 0._wp 
     393                  END DO 
     394               END DO 
     395            ENDIF 
     396            ! 
     397         END DO 
     398         ! 
     399      ENDIF  
     400#endif 
     401      ! 
     402      IF( ln_timing )   CALL timing_stop('wzv_RK3') 
     403      ! 
     404   END SUBROUTINE wzv_RK3 
    261405 
    262406 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/eosbn2.F90

    r14131 r14547  
    5656   !                  !! * Interface 
    5757   INTERFACE eos 
    58       MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d 
     58      MODULE PROCEDURE eos_insitu_New, eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d 
    5959   END INTERFACE 
    6060   ! 
     
    188188   !!---------------------------------------------------------------------- 
    189189CONTAINS 
     190 
     191   SUBROUTINE eos_insitu_New( pts, Knn, prd ) 
     192      !!---------------------------------------------------------------------- 
     193      !!                   ***  ROUTINE eos_insitu  *** 
     194      !! 
     195      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
     196      !!       potential temperature and salinity using an equation of state 
     197      !!       selected in the nameos namelist 
     198      !! 
     199      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 
     200      !!         with   prd    in situ density anomaly      no units 
     201      !!                t      TEOS10: CT or EOS80: PT      Celsius 
     202      !!                s      TEOS10: SA or EOS80: SP      TEOS10: g/kg or EOS80: psu 
     203      !!                z      depth                        meters 
     204      !!                rho    in situ density              kg/m^3 
     205      !!                rho0   reference density            kg/m^3 
     206      !! 
     207      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 
     208      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg 
     209      !! 
     210      !!     ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z). 
     211      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu 
     212      !! 
     213      !!     ln_seos : simplified equation of state 
     214      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 
     215      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0 
     216      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 
     217      !!              Vallis like equation: use default values of coefficients 
     218      !! 
     219      !! ** Action  :   compute prd , the in situ density (no units) 
     220      !! 
     221      !! References :   Roquet et al, Ocean Modelling, in preparation (2014) 
     222      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     223      !!                TEOS-10 Manual, 2010 
     224      !!---------------------------------------------------------------------- 
     225      REAL(wp), DIMENSION(:,:,:,:,:), INTENT(in   ) ::   pts   ! T-S 
     226      INTEGER                     , INTENT(in   ) ::   Knn   ! time-level 
     227      REAL(wp), DIMENSION(:,:,:  ), INTENT(  out) ::   prd   ! in situ density 
     228      ! 
     229      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     230      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars 
     231      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     232      !!---------------------------------------------------------------------- 
     233      ! 
     234      IF( ln_timing )   CALL timing_start('eos-insitu') 
     235      ! 
     236      SELECT CASE( neos ) 
     237      ! 
     238      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     239         ! 
     240         DO_3D(nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     241            ! 
     242            zh  = gdept(ji,jj,jk,Knn) * r1_Z0                                 ! depth 
     243            zt  = pts (ji,jj,jk,jp_tem,Knn) * r1_T0                           ! temperature 
     244            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal,Knn) + rdeltaS ) * r1_S0 )   ! square root salinity 
     245            ztm = tmask(ji,jj,jk)                                             ! tmask 
     246            ! 
     247            zn3 = EOS013*zt   & 
     248               &   + EOS103*zs+EOS003 
     249               ! 
     250            zn2 = (EOS022*zt   & 
     251               &   + EOS112*zs+EOS012)*zt   & 
     252               &   + (EOS202*zs+EOS102)*zs+EOS002 
     253               ! 
     254            zn1 = (((EOS041*zt   & 
     255               &   + EOS131*zs+EOS031)*zt   & 
     256               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     257               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     258               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     259               ! 
     260            zn0 = (((((EOS060*zt   & 
     261               &   + EOS150*zs+EOS050)*zt   & 
     262               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     263               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     264               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     265               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     266               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     267               ! 
     268            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     269            ! 
     270            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked) 
     271            ! 
     272         END_3D 
     273         ! 
     274      CASE( np_seos )                !==  simplified EOS  ==! 
     275         ! 
     276         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     277            zt  = pts  (ji,jj,jk,jp_tem,Knn) - 10._wp 
     278            zs  = pts  (ji,jj,jk,jp_sal,Knn) - 35._wp 
     279            zh  = gdept(ji,jj,jk,Knn) 
     280            ztm = tmask(ji,jj,jk) 
     281            ! 
     282            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     283               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     284               &  - rn_nu * zt * zs 
     285               ! 
     286            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
     287         END_3D 
     288         ! 
     289      END SELECT 
     290      ! 
     291      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk ) 
     292      ! 
     293      IF( ln_timing )   CALL timing_stop('eos-insitu') 
     294      ! 
     295   END SUBROUTINE eos_insitu_New 
     296 
    190297 
    191298   SUBROUTINE eos_insitu( pts, prd, pdep ) 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/tranpc.F90

    r14215 r14547  
    7171      LOGICAL  ::   l_bottom_reached, l_column_treated 
    7272      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 
    73       REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt 
     73      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw 
    7474      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp             ! acceptance criteria for neutrality (N2==0) 
    7575      REAL(wp), DIMENSION(    jpk     )   ::   zvn2             ! vertical profile of N2 at 1 given point... 
     
    311311         ! 
    312312         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic 
    313             z1_rDt = 1._wp / (2._wp * rn_Dt) 
    314             ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt 
    315             ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt 
     313            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * r1_Dt 
     314            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * r1_Dt 
    316315            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt ) 
    317316            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds ) 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/nemogcm.F90

    r14239 r14547  
    6767#endif 
    6868#if defined key_qco   ||   defined key_linssh 
     69# if defined key_RK3 
     70   USE stprk3 
     71# else 
    6972   USE stpmlf         ! NEMO time-stepping               (stp_MLF   routine) 
     73# endif 
    7074#else 
    7175   USE step           ! NEMO time-stepping                 (stp     routine) 
     
    163167         ! 
    164168#  if defined key_qco   ||   defined key_linssh 
     169#   if defined key_RK3 
     170         CALL stp_RK3 
     171#   else 
    165172         CALL stp_MLF 
     173#   endif 
    166174#  else 
    167175         CALL stp 
     
    184192            ! 
    185193#  if defined key_qco   ||   defined key_linssh 
     194#   if defined key_RK3 
     195            CALL stp_RK3( istp ) 
     196#   else 
    186197            CALL stp_MLF( istp ) 
     198#   endif 
    187199#  else 
    188200            CALL stp    ( istp ) 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90

    r14381 r14547  
    6262#  include "do_loop_substitute.h90" 
    6363#  include "domzgr_substitute.h90" 
    64 #  include "do_loop_substitute.h90" 
    6564   !!---------------------------------------------------------------------- 
    6665   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
Note: See TracChangeset for help on using the changeset viewer.