Changeset 13603


Ignore:
Timestamp:
2020-10-14T18:04:16+02:00 (3 months ago)
Author:
techene
Message:

#2385 adapted for SWE : some cleaning, bug corrected in the time stepping (more details in the ticket)

File:
1 moved

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpMLF.F90

    r13602 r13603  
    1 MODULE stpLF 
     1MODULE stpMLF 
    22   !!====================================================================== 
    3    !!                       ***  MODULE step  *** 
     3   !!                       ***  MODULE stpMLF  *** 
    44   !! Time-stepping   : manager of the shallow water equation time stepping 
     5   !!                   Modified Leap Frog scheme 
    56   !!====================================================================== 
    67   !! History :  NEMO !  2020-03  (A. Nasser, G. Madec)  Original code from  4.0.2 
    7    !!---------------------------------------------------------------------- 
    8  
    9    !!---------------------------------------------------------------------- 
    10    !!   stp             : Shallow Water time-stepping 
    11    !!---------------------------------------------------------------------- 
    12    USE step_oce         ! time stepping definition modules 
    13    USE phycst           ! physical constants 
    14    USE usrdef_nam 
     8   !!             -   !  2020-10  (S. Techene, G. Madec)  cleanning 
     9   !!---------------------------------------------------------------------- 
     10 
     11   !!---------------------------------------------------------------------- 
     12   !!   stp_MLF       : MLF Shallow Water Eq. time-stepping 
     13   !!---------------------------------------------------------------------- 
     14   USE stp_oce        ! modules used in nemo_init and stp_MLF 
    1515   ! 
    16    USE iom              ! xIOs server  
    17    USE domqco 
    18  
     16   USE domqco         ! quasi-eulerian coordinate 
     17   USE phycst         ! physical constants 
     18   USE usrdef_nam     ! user defined namelist parameters 
     19!!st   USE usrdef_sbc     ! user defined surface boundary cond 
     20    
    1921   IMPLICIT NONE 
    2022   PRIVATE 
    2123 
    22    PUBLIC   stp_LF   ! called by nemogcm.F90 
     24   PUBLIC   stp_MLF   ! called by nemogcm.F90 
    2325    
    24    !!---------------------------------------------------------------------- 
    25    !! time level indices 
    26    !!---------------------------------------------------------------------- 
    27    INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !! used by nemo_init 
     26   !                                          !**  time level indices  **! 
     27   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !: used by nemo_init 
    2828       
    2929   !! * Substitutions 
     
    3737CONTAINS 
    3838 
    39 #if defined key_agrif 
    40    RECURSIVE SUBROUTINE stp_LF( ) 
    41       INTEGER             ::   kstp   ! ocean time-step index 
    42 #else 
    43    SUBROUTINE stp_LF( kstp ) 
    44       INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
    45 #endif 
     39   SUBROUTINE stp_MLF( kstp ) 
    4640      !!---------------------------------------------------------------------- 
    47       !!                     ***  ROUTINE stp  *** 
     41      !!                     ***  ROUTINE stp_MLF  *** 
    4842      !! 
    4943      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.) 
     
    5650      !!              -6- Outputs and diagnostics 
    5751      !!---------------------------------------------------------------------- 
    58       INTEGER ::   ji, jj, jk   ! dummy loop indice 
    59       INTEGER ::   indic        ! error indicator if < 0 
    60 !!gm kcall can be removed, I guess 
    61       INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt) 
    62       REAL(wp)::   z1_2rho0     ! local scalars 
    63        
    64       REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars 
    65       REAL(wp) ::   zve3a, zve3n, zve3b    !   -      - 
    66       REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva 
     52      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
     53      ! 
     54      INTEGER ::   ji, jj, jk             ! dummy loop indices 
     55      INTEGER ::   indic                  ! error indicator if < 0 
     56      REAL(wp)::   z1_2rho0               ! local scalars 
     57      REAL(wp)::   zrhs_u, zue3a, zue3n, zue3b, zua   ! local scalars 
     58      REAL(wp)::   zrhs_v, zve3a, zve3n, zve3b, zva   !   -      - 
    6759      !! --------------------------------------------------------------------- 
    68 #if defined key_agrif 
    69       kstp = nit000 + Agrif_Nb_Step() 
    70       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices 
    71       IF( lk_agrif_debug ) THEN 
    72          IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---' 
    73          IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint() 
    74       ENDIF 
    75       IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE. 
    76 # if defined key_iomput 
    77       IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context ) 
    78 # endif 
    79 #endif 
    80       ! 
    81       IF( ln_timing )   CALL timing_start('stp_LF') 
     60      ! 
     61      IF( ln_timing )   CALL timing_start('stp_MLF') 
    8262      ! 
    8363      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    8565      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    8666      ! 
    87       IF( l_1st_euler ) THEN     ! start or restart with Euler 1st time-step 
    88          rDt   =  rn_Dt    
     67      IF( l_1st_euler ) THEN   ! start or restart with Euler 1st time-step 
     68         rDt   = rn_Dt    
    8969         r1_Dt = 1._wp / rDt 
    9070      ENDIF 
    9171 
    92       IF ( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero 
     72      IF( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero 
    9373 
    9474      ! 
     
    9777      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    9878                             indic = 0                ! reset to no error condition 
    99                               
     79 
    10080      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS) 
    10181                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom) 
    102          IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis 
    10382                             CALL iom_init_closedef 
    104          IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid 
    10583      ENDIF 
    10684      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init) 
    10785                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp 
    108       IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp 
    109  
    110       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    111       ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice) 
    112       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    113       IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential 
    114       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)  
    115       IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries 
    116                          CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition 
    117  
    118       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    119       ! Ocean physics update 
    120       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    121       !  LATERAL  PHYSICS 
    122       !                                                                        ! eddy diffusivity coeff. 
     86 
     87      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     88      ! Update external forcing   (SWE: surface boundary condition only) 
     89      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     90 
     91                             CALL sbc     ( kstp, Nbb, Nnn )                   ! Sea Boundary Condition 
     92 
     93      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     94      ! Ocean physics update   (SWE: eddy viscosity only) 
     95      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     96 
    12397      IF( l_ldfdyn_time )   CALL ldf_dyn( kstp, Nbb )                          ! eddy viscosity coeff.  
    12498 
    12599      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    126       !  Ocean dynamics : hdiv, ssh, e3, u, v, w 
    127       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    128  
    129                             CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor) 
     100      !  RHS of horizontal velocity 
     101      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     102 
    130103                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero 
    131104                         vv(:,:,:,Nrhs) = 0._wp 
    132105 
    133       IF( ln_bdy     )      CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends 
    134  
    135 #if defined key_agrif 
    136       IF(.NOT. Agrif_Root())  &  
    137                &            CALL Agrif_Sponge_dyn        ! momentum sponge 
    138 #endif 
    139                             CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
    140   
    141                             CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    142   
    143                             CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
    144  
    145       IF( .NOT.ln_linssh )  CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit 
    146       !IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt_st( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors  
    147 !!an - calcul du gradient de pression horizontal (explicit) 
     106                         CALL dyn_adv( kstp, Nbb, Nnn, uu, vv, Nrhs )   ! advection (VF or FF)  ==> RHS 
     107                         CALL dyn_vor( kstp,      Nnn, uu, vv, Nrhs )   ! vorticity             ==> RHS 
     108                         CALL dyn_ldf( kstp, Nbb, Nnn, uu, vv, Nrhs )   ! lateral mixing                     ==> RHS 
     109 
     110      z1_2rho0 = 0.5_wp * r1_rho0 
     111      ! 
    148112      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    149          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 
    150          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 
     113         !                                                              ! horizontal pressure gradient 
     114         zrhs_u =        - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 
     115         zrhs_v =        - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 
     116         !                                                              ! wind stress and layer friction 
     117         zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   & 
     118            &                                  - rn_rfr * uu(ji,jj,jk,Nbb) 
     119         zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   & 
     120            &                                  - rn_rfr * vv(ji,jj,jk,Nbb) 
     121         !                                                              ! ==> RHS 
     122         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u 
     123         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v 
    151124      END_3D 
    152       ! 
    153       ! add wind stress forcing and layer linear friction to the RHS  
    154       z1_2rho0 = 0.5_wp * r1_rho0 
    155       DO_3D( 0, 0, 0, 0,1,jpkm1) 
    156          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   & 
    157             &                                  - rn_rfr * uu(ji,jj,jk,Nbb) 
    158          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   & 
    159             &                                  - rn_rfr * vv(ji,jj,jk,Nbb) 
    160       END_3D    
    161 !!an          
    162    
    163       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    164       ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3  
    165       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    166 !!st ssh_atf : add ssh filtering up there 
    167       IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps 
    168          !                                                  ! filtering "now" field 
     125 
     126      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     127      ! Time stepping of ssh Eq. (and update r3_Naa) 
     128      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     129      !        ! Leap Frog time stepping ==> ssh_Naa and r3_Naa 
     130                         CALL ssh_nxt( kstp, Nbb, Nnn, ssh   , Naa )    ! after ssh  
     131      !                                                                 ! after ssh/h_0 ratio explicit 
     132                         CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) ) 
     133      !        ! Asselin filter ==> ssh_Nnn filtered 
     134      IF ( .NOT.( l_1st_euler ) ) THEN   ! Time filtering of now ssh 
    169135         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2._wp * ssh(:,:,Nnn) + ssh(:,:,Naa) ) 
    170136      ENDIF 
    171 !!st ssh_atf 
    172            
    173 !! what about  IF( .NOT.ln_linssh )  ? 
    174 !!an futur module dyn_nxt (a la place de dyn_atf) 
    175        
     137 
     138      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     139      ! Time stepping of dynamics (u,v) 
     140      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     141 
    176142      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity 
    177          IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter) 
    178             DO_3D( 0, 0, 0, 0,1,jpkm1) 
     143         IF( l_1st_euler ) THEN          ! Euler time stepping (no Asselin filter) 
     144            DO_3D( 0,0, 0,0, 1,jpkm1) 
    179145               uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
    180146               vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
    181147             END_3D           
    182          ELSE                             ! Leap Frog time stepping + Asselin filter          
    183             DO_3D( 1, 1, 1, 1,1,jpkm1) 
     148         ELSE                            ! Leap Frog time stepping + Asselin filter          
     149            DO_3D( 0,0, 0,0, 1,jpkm1) 
    184150               zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
    185151               zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     
    191157               vv(ji,jj,jk,Naa) = zva 
    192158            END_3D 
    193             CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn) )   ! "now" ssh/h_0 ratio from filtrered ssh 
    194 #if ! defined key_qco 
    195             DO jk = 1, jpkm1 
    196                e3t(:,:,jk,Nnn) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Nnn) ) 
    197                e3u(:,:,jk,Nnn) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Nnn) ) 
    198                e3v(:,:,jk,Nnn) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Nnn) ) 
    199             END DO 
    200 #endif 
     159            !                            ! Update r3_Nnn 
     160            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn) )   ! now ssh/h_0 ratio from filtrered ssh 
    201161         ENDIF 
    202162         ! 
    203163      ELSE                          ! flux form : applied on thickness weighted velocity 
    204          IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter) 
    205             DO_3D( 0, 0, 0, 0,1,jpkm1) 
     164         IF( l_1st_euler ) THEN          ! Euler time stepping (no Asselin filter) 
     165            DO_3D( 0,0, 0,0, 1,jpkm1) 
    206166               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 
    207167               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 
    208                !                                                ! LF time stepping 
    209                zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
    210                zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     168               !                                                ! Euler time stepping 
     169               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     170               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
    211171               ! 
    212172               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)     
     
    214174            END_3D 
    215175         ELSE                             ! Leap Frog time stepping + Asselin filter 
    216             CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:) )   ! "now" ssh/h_0 ratio from filtrered ssh 
    217             DO_3D( 1, 1, 1, 1,1,jpkm1) 
    218                zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn) 
    219                zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn) 
    220                zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 
    221                zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 
     176            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:) )   ! now ssh/h_0 ratio from filtrered ssh 
     177            DO_3D( 0,0, 0,0, 1,jpkm1) 
     178               zue3n = ( 1._wp + r3u(ji,jj,Nnn) ) * uu(ji,jj,jk,Nnn) 
     179               zve3n = ( 1._wp + r3v(ji,jj,Nnn) ) * vv(ji,jj,jk,Nnn) 
     180               zue3b = ( 1._wp + r3u(ji,jj,Nbb) ) * uu(ji,jj,jk,Nbb) 
     181               zve3b = ( 1._wp + r3v(ji,jj,Nbb) ) * vv(ji,jj,jk,Nbb) 
    222182               !                                                ! LF time stepping 
    223                zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
    224                zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
    225                !                                                ! Asselin time filter on e3u/v/t 
    226                ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn)  + e3u(ji,jj,jk,Naa) )  
    227                ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn)  + e3v(ji,jj,jk,Naa) ) 
    228                ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn)  + e3t(ji,jj,jk,Naa) )  
     183               zue3a = zue3b + rDt * ( 1._wp + r3u(ji,jj,Nnn) ) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     184               zve3a = zve3b + rDt * ( 1._wp + r3v(ji,jj,Nnn) ) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
    229185               !                                                ! Asselin time filter on u,v (Nnn) 
    230                uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / (e3u_0(ji,jj,jk) * ( 1._wp + r3u_f(ji,jj) )) 
    231                vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / (e3v_0(ji,jj,jk) * ( 1._wp + r3v_f(ji,jj) )) 
     186               uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ( 1._wp + r3u_f(ji,jj) ) 
     187               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ( 1._wp + r3v_f(ji,jj) ) 
    232188               ! 
    233                uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)     
    234                vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 
     189               uu(ji,jj,jk,Naa) = zue3a / ( 1._wp + r3u(ji,jj,Naa) )     
     190               vv(ji,jj,jk,Naa) = zve3a / ( 1._wp + r3v(ji,jj,Naa) ) 
    235191            END_3D 
     192            !                             ! Update r3_Nnn with time filtered values 
    236193            r3t(:,:,Nnn) = r3t_f(:,:) 
    237194            r3u(:,:,Nnn) = r3u_f(:,:) 
    238195            r3v(:,:,Nnn) = r3v_f(:,:) 
    239 #if ! defined key_qco 
    240             DO jk = 1, jpkm1 
    241                e3t(:,:,jk,Nnn) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Nnn) ) 
    242                e3u(:,:,jk,Nnn) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Nnn) ) 
    243                e3v(:,:,jk,Nnn) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Nnn) ) 
    244             END DO 
    245 #endif          
    246196         ENDIF 
    247197      ENDIF 
    248198 
    249  
    250       CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries 
    251          &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )      
    252  
    253 !!an          
     199      CALL lbc_lnk_multi( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries 
     200         &                           uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )      
    254201 
    255202      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    256203      ! Set boundary conditions, time filter and swap time levels 
    257204      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    258  
    259 !!an TO BE ADDED : dyn_nxt 
    260 !!                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors 
    261 !!an TO BE ADDED : a simplifier 
    262 !                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height 
    263 !!st copied above  
    264 !!      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps 
    265 !!         !                                                  ! filtering "now" field 
    266 !!         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) ) 
    267 !!      ENDIF 
    268 !!st       
    269 !!an  
    270  
    271205 
    272206      ! Swap time levels 
     
    275209      Nnn = Naa 
    276210      Naa = Nrhs 
    277       ! 
    278 !                         CALL dom_vvl_sf_update_st( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors 
     211 
    279212      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    280213      ! diagnostics and outputs 
    281214      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    282       IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats 
    283       IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics 
    284      
    285                          CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs 
    286  
    287       ! 
    288       IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file 
    289            
    290  
    291 #if defined key_agrif 
    292       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    293       ! AGRIF 
    294       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<       
    295                          Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices 
    296                          CALL Agrif_Integrate_ChildGrids( stp_LF )       ! allows to finish all the Child Grids before updating 
    297  
    298                          IF( Agrif_NbStepint() == 0 ) THEN 
    299                             CALL Agrif_update_all( )                  ! Update all components 
    300                          ENDIF 
    301 #endif 
    302       IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
     215       
     216      IF( ln_diacfl  )   CALL dia_cfl  ( kstp,      Nnn )      ! Courant number diagnostics 
     217                         CALL dia_wri  ( kstp,      Nnn )      ! ocean model: outputs 
     218      ! 
     219      IF( lrst_oce   )   CALL rst_write( kstp, Nbb, Nnn )      ! write output ocean restart file 
    303220 
    304221      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    305222      ! Control 
    306223      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    307                          CALL stp_ctl      ( kstp, Nbb, Nnn, indic ) 
    308                 
    309                           
     224                         CALL stp_ctl      ( kstp, Nnn ) 
     225 
    310226      IF( kstp == nit000 ) THEN                          ! 1st time step only 
    311227                                        CALL iom_close( numror )   ! close input  ocean restart file 
     
    316232      ! 
    317233#if defined key_iomput 
    318       IF( kstp == nitend .OR. indic < 0 ) THEN  
    319                       CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF 
    320                       IF(lrxios) CALL iom_context_finalize(      crxios_context          ) 
     234      IF( kstp == nitend .OR. indic < 0 ) THEN 
     235!!st : cxios_context needed ? because opened earlier ???          
     236         CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF 
     237!!st : crxios_context not needed associated with coarsening !  
     238         IF(lrxios) CALL iom_context_finalize( crxios_context ) 
    321239      ENDIF 
    322240#endif 
    323241      ! 
    324242      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep 
    325          rDt = 2._wp * rn_Dt    
     243         rDt   = 2._wp * rn_Dt    
    326244         r1_Dt = 1._wp / rDt 
    327245         l_1st_euler = .FALSE.       
    328246      ENDIF 
    329247      ! 
    330       IF( ln_timing )   CALL timing_stop('stp') 
    331       ! 
    332    END SUBROUTINE stp_LF 
    333    ! 
     248      IF( ln_timing )   CALL timing_stop('stp_MLF') 
     249      ! 
     250   END SUBROUTINE stp_MLF 
     251 
    334252   !!====================================================================== 
    335 END MODULE stpLF 
     253END MODULE stpMLF 
Note: See TracChangeset for help on using the changeset viewer.