Changeset 12581


Ignore:
Timestamp:
2020-03-20T23:26:56+01:00 (3 weeks ago)
Author:
techene
Message:

OCE/DOM/domqe.F90: make dom_qe_r3c public, OCE/DYN/dynatfLF.F90: duplicate dynatf and replace dom_qe_interpol calls by the ssh scaling method to compute and update e3t/u/v at time Kmm, OCE/TRA/traatfLF.F90: duplicate traatf and replace dom_qe_interpol by the ssh scaling method to compute internal ze3t, OCE/steplf.F90: change the order of atf routines ssh_atf is called first then dom_qe_r3c to computed filtered ssh to h ratios at T-, U-, V-points finally tra_atf_lf and dyn_atf_lf

Location:
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE
Files:
4 edited
2 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90

    r12579 r12581  
    4545   PUBLIC  dom_qe_sf_update  ! called by step.F90 
    4646   PUBLIC  dom_qe_interpol   ! called by dynnxt.F90 
     47   PUBLIC  dom_qe_r3c        ! called by steplf.F90 
    4748 
    4849   !                                                      !!* Namelist nam_vvl 
  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90

    r12377 r12581  
    1313   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond. 
    1414   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
    15    !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.  
     15   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines. 
    1616   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option 
    1717   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module 
     
    2222   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering. 
    2323   !!------------------------------------------------------------------------- 
    24    
     24 
    2525   !!---------------------------------------------------------------------------------------------- 
    2626   !!   dyn_atf       : apply Asselin time filtering to "now" velocities and vertical scale factors 
     
    4242   USE trdken         ! trend manager: kinetic energy 
    4343   USE isf_oce   , ONLY: ln_isf     ! ice shelf 
    44    USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine  
     44   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 
    4545   ! 
    4646   USE in_out_manager ! I/O manager 
     
    6363   !!---------------------------------------------------------------------- 
    6464   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    65    !! $Id$  
     65   !! $Id$ 
    6666   !! Software governed by the CeCILL license (see ./LICENSE) 
    6767   !!---------------------------------------------------------------------- 
     
    7171      !!---------------------------------------------------------------------- 
    7272      !!                  ***  ROUTINE dyn_atf  *** 
    73       !!                    
    74       !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary  
     73      !! 
     74      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary 
    7575      !!             condition on the after velocity and apply the Asselin time 
    7676      !!             filter to the now fields. 
     
    7979      !!             estimate (ln_dynspg_ts=T) 
    8080      !! 
    81       !!              * Apply lateral boundary conditions on after velocity  
     81      !!              * Apply lateral boundary conditions on after velocity 
    8282      !!             at the local domain boundaries through lbc_lnk call, 
    8383      !!             at the one-way open boundaries (ln_bdy=T), 
     
    8686      !!              * Apply the Asselin time filter to the now fields 
    8787      !!             arrays to start the next time step: 
    88       !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))  
     88      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 
    8989      !!                                    + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ] 
    9090      !!             Note that with flux form advection and non linear free surface, 
     
    9292      !!             As a result, dyn_atf MUST be called after tra_atf. 
    9393      !! 
    94       !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity  
     94      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity 
    9595      !!---------------------------------------------------------------------- 
    9696      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index 
     
    103103      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      - 
    104104      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zwfld 
    105       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva  
     105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva 
    106106      !!---------------------------------------------------------------------- 
    107107      ! 
     
    131131         ! 
    132132         IF( .NOT.ln_bt_fw ) THEN 
    133             ! Remove advective velocity from "now velocities"  
    134             ! prior to asselin filtering      
    135             ! In the forward case, this is done below after asselin filtering    
    136             ! so that asselin contribution is removed at the same time  
     133            ! Remove advective velocity from "now velocities" 
     134            ! prior to asselin filtering 
     135            ! In the forward case, this is done below after asselin filtering 
     136            ! so that asselin contribution is removed at the same time 
    137137            DO jk = 1, jpkm1 
    138138               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 
    139139               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 
    140             END DO   
     140            END DO 
    141141         ENDIF 
    142142      ENDIF 
    143143 
    144144      ! Update after velocity on domain lateral boundaries 
    145       ! --------------------------------------------------       
     145      ! -------------------------------------------------- 
    146146# if defined key_agrif 
    147147      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries 
     
    157157      ! 
    158158      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics 
    159          z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step  
     159         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step 
    160160         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt 
    161161         ! 
     
    177177      ! Time filter and swap of dynamics arrays 
    178178      ! ------------------------------------------ 
    179           
    180       IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter  
     179 
     180      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter 
    181181         !                                ! =============! 
    182182         IF( ln_linssh ) THEN             ! Fixed volume ! 
     
    202202            DO jk = 1, jpkm1 
    203203               ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) & 
    204                               &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) )  
     204                              &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 
    205205            END DO 
    206206            ! 
     
    239239                  pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
    240240               END_3D 
    241                pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1)   
     241               pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 
    242242               pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 
    243243               ! 
     
    250250         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    251251            ! Revert filtered "now" velocities to time split estimate 
    252             ! Doing it here also means that asselin filter contribution is removed   
     252            ! Doing it here also means that asselin filter contribution is removed 
    253253            zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 
    254             zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)     
     254            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 
    255255            DO jk = 2, jpkm1 
    256256               zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
    257                zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)     
     257               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 
    258258            END DO 
    259259            DO jk = 1, jpkm1 
     
    307307      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   & 
    308308         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask ) 
    309       !  
     309      ! 
    310310      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve ) 
    311311      IF( l_trddyn     )   DEALLOCATE( zua, zva ) 
  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90

    r12481 r12581  
    1 MODULE dynatf 
     1MODULE dynatfLF 
    22   !!========================================================================= 
    3    !!                       ***  MODULE  dynatf  *** 
     3   !!                       ***  MODULE  dynatfLF  *** 
    44   !! Ocean dynamics: time filtering 
    55   !!========================================================================= 
     
    1313   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond. 
    1414   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
    15    !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.  
     15   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines. 
    1616   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option 
    1717   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module 
     
    2020   !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends 
    2121   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification 
    22    !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering. 
     22   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatfLF.F90. Now just does time filtering. 
    2323   !!------------------------------------------------------------------------- 
    24    
     24 
    2525   !!---------------------------------------------------------------------------------------------- 
    26    !!   dyn_atf       : apply Asselin time filtering to "now" velocities and vertical scale factors 
     26   !!   dyn_atfLF       : apply Asselin time filtering to "now" velocities and vertical scale factors 
    2727   !!---------------------------------------------------------------------------------------------- 
    2828   USE oce            ! ocean dynamics and tracers 
     
    4242   USE trdken         ! trend manager: kinetic energy 
    4343   USE isf_oce   , ONLY: ln_isf     ! ice shelf 
    44    USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine  
     44   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 
    4545   ! 
    4646   USE in_out_manager ! I/O manager 
     
    5757   PRIVATE 
    5858 
    59    PUBLIC    dyn_atf   ! routine called by step.F90 
     59   PUBLIC    dyn_atf_lf   ! routine called by step.F90 
    6060 
    6161   !! * Substitutions 
     
    6363   !!---------------------------------------------------------------------- 
    6464   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    65    !! $Id$  
     65   !! $Id$ 
    6666   !! Software governed by the CeCILL license (see ./LICENSE) 
    6767   !!---------------------------------------------------------------------- 
    6868CONTAINS 
    6969 
    70    SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 
     70   SUBROUTINE dyn_atf_lf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 
    7171      !!---------------------------------------------------------------------- 
    72       !!                  ***  ROUTINE dyn_atf  *** 
    73       !!                    
    74       !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary  
     72      !!                  ***  ROUTINE dyn_atf_lf  *** 
     73      !! 
     74      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary 
    7575      !!             condition on the after velocity and apply the Asselin time 
    7676      !!             filter to the now fields. 
     
    7979      !!             estimate (ln_dynspg_ts=T) 
    8080      !! 
    81       !!              * Apply lateral boundary conditions on after velocity  
     81      !!              * Apply lateral boundary conditions on after velocity 
    8282      !!             at the local domain boundaries through lbc_lnk call, 
    8383      !!             at the one-way open boundaries (ln_bdy=T), 
     
    8686      !!              * Apply the Asselin time filter to the now fields 
    8787      !!             arrays to start the next time step: 
    88       !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))  
     88      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 
    8989      !!                                    + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ] 
    9090      !!             Note that with flux form advection and non linear free surface, 
    9191      !!             the time filter is applied on thickness weighted velocity. 
    92       !!             As a result, dyn_atf MUST be called after tra_atf. 
    93       !! 
    94       !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity  
     92      !!             As a result, dyn_atf_lf MUST be called after tra_atf. 
     93      !! 
     94      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity 
    9595      !!---------------------------------------------------------------------- 
    9696      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index 
     
    102102      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars 
    103103      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      - 
    104       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zwfld 
    105       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva  
     104      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
     105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zua, zva 
    106106      !!---------------------------------------------------------------------- 
    107107      ! 
    108       IF( ln_timing    )   CALL timing_start('dyn_atf') 
     108      IF( ln_timing    )   CALL timing_start('dyn_atf_lf') 
    109109      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     ) 
    110110      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) ) 
     
    112112      IF( kt == nit000 ) THEN 
    113113         IF(lwp) WRITE(numout,*) 
    114          IF(lwp) WRITE(numout,*) 'dyn_atf : Asselin time filtering' 
     114         IF(lwp) WRITE(numout,*) 'dyn_atf_lf : Asselin time filtering' 
    115115         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    116116      ENDIF 
     
    131131         ! 
    132132         IF( .NOT.ln_bt_fw ) THEN 
    133             ! Remove advective velocity from "now velocities"  
    134             ! prior to asselin filtering      
    135             ! In the forward case, this is done below after asselin filtering    
    136             ! so that asselin contribution is removed at the same time  
     133            ! Remove advective velocity from "now velocities" 
     134            ! prior to asselin filtering 
     135            ! In the forward case, this is done below after asselin filtering 
     136            ! so that asselin contribution is removed at the same time 
    137137            DO jk = 1, jpkm1 
    138138               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 
    139139               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 
    140             END DO   
     140            END DO 
    141141         ENDIF 
    142142      ENDIF 
    143143 
    144144      ! Update after velocity on domain lateral boundaries 
    145       ! --------------------------------------------------       
     145      ! -------------------------------------------------- 
    146146# if defined key_agrif 
    147147      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries 
    148148# endif 
    149149      ! 
    150       CALL lbc_lnk_multi( 'dynatf', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. )     !* local domain boundaries 
     150      CALL lbc_lnk_multi( 'dynatfLF', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. )     !* local domain boundaries 
    151151      ! 
    152152      !                                !* BDY open boundaries 
     
    157157      ! 
    158158      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics 
    159          z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step  
     159         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step 
    160160         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt 
    161161         ! 
     
    177177      ! Time filter and swap of dynamics arrays 
    178178      ! ------------------------------------------ 
    179           
    180       IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter  
     179 
     180      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter 
    181181         !                                ! =============! 
    182182         IF( ln_linssh ) THEN             ! Fixed volume ! 
     
    191191            ! Time-filtered scale factor at t-points 
    192192            ! ---------------------------------------------------- 
    193             ALLOCATE( ze3t_f(jpi,jpj,jpk), zwfld(jpi,jpj) ) 
    194             DO jk = 1, jpkm1 
    195                ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) ) 
    196             END DO 
    197             ! Add volume filter correction: compatibility with tracer advection scheme 
    198             ! => time filter + conservation correction 
    199             zcoef = atfp * rdt * r1_rau0 
    200             zwfld(:,:) = emp_b(:,:) - emp(:,:) 
    201             IF ( ln_rnf ) zwfld(:,:) =  zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) ) 
    202             DO jk = 1, jpkm1 
    203                ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) & 
    204                               &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) )  
     193            DO jk = 1, jpk                                          ! filtered scale factor at T-points 
     194               pe3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t_f(:,:) * tmask(:,:,jk) ) 
    205195            END DO 
    206196            ! 
    207             ! ice shelf melting (deal separately as it can be in depth) 
    208             ! PM: we could probably define a generic subroutine to do the in depth correction 
    209             !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...) 
    210             !     ...(kt, coef, ktop, kbot, hz, fwf_b, fwf) 
    211             IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, atfp * rdt ) 
    212             ! 
    213             pe3t(:,:,1:jpkm1,Kmm) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points 
    214197            ! 
    215198            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 
    216199               ! Before filtered scale factor at (u/v)-points 
    217                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 
    218                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 
     200               DO jk = 1, jpk 
     201                  pe3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u_f(:,:) * umask(:,:,jk) ) 
     202                  pe3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v_f(:,:) * vmask(:,:,jk) ) 
     203               END DO 
     204               ! 
    219205               DO_3D_11_11( 1, jpkm1 ) 
    220206                  puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 
     
    224210            ELSE                          ! Asselin filter applied on thickness weighted velocity 
    225211               ! 
    226                ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 
    227                ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 
    228                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' ) 
    229                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' ) 
    230212               DO_3D_11_11( 1, jpkm1 ) 
    231213                  zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) 
     
    235217                  zue3b = pe3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb) 
    236218                  zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb) 
     219                  !                                                 ! filtered scale factor at U-,V-points 
     220                  pe3u(ji,jj,jk,Kmm) = e3u_0(ji,jj,jk) * ( 1._wp + r3u_f(ji,jj) * umask(ji,jj,jk) ) 
     221                  pe3v(ji,jj,jk,Kmm) = e3v_0(ji,jj,jk) * ( 1._wp + r3v_f(ji,jj) * vmask(ji,jj,jk) ) 
    237222                  ! 
    238                   puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
    239                   pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
     223                  puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / pe3u(ji,jj,jk,Kmm) !!st ze3u_f(ji,jj,jk) 
     224                  pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / pe3v(ji,jj,jk,Kmm) !!st ze3v_f(ji,jj,jk) 
    240225               END_3D 
    241                pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1)   
    242                pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 
    243                ! 
    244                DEALLOCATE( ze3u_f , ze3v_f ) 
     226               ! 
    245227            ENDIF 
    246228            ! 
    247             DEALLOCATE( ze3t_f, zwfld ) 
    248229         ENDIF 
    249230         ! 
    250231         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    251232            ! Revert filtered "now" velocities to time split estimate 
    252             ! Doing it here also means that asselin filter contribution is removed   
     233            ! Doing it here also means that asselin filter contribution is removed 
    253234            zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 
    254             zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)     
     235            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 
    255236            DO jk = 2, jpkm1 
    256237               zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
    257                zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)     
     238               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 
    258239            END DO 
    259240            DO jk = 1, jpkm1 
     
    307288      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   & 
    308289         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask ) 
    309       !  
     290      ! 
    310291      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve ) 
    311292      IF( l_trddyn     )   DEALLOCATE( zua, zva ) 
    312       IF( ln_timing    )   CALL timing_stop('dyn_atf') 
    313       ! 
    314    END SUBROUTINE dyn_atf 
     293      IF( ln_timing    )   CALL timing_stop('dyn_atf_lf') 
     294      ! 
     295   END SUBROUTINE dyn_atf_lf 
    315296 
    316297   !!========================================================================= 
    317 END MODULE dynatf 
     298END MODULE dynatfLF 
  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90

    r12377 r12581  
    2626   !!---------------------------------------------------------------------- 
    2727   USE oce             ! ocean dynamics and tracers variables 
    28    USE dom_oce         ! ocean space and time domain variables  
     28   USE dom_oce         ! ocean space and time domain variables 
    2929   USE sbc_oce         ! surface boundary condition: ocean 
    3030   USE sbcrnf          ! river runoffs 
     
    3333   USE domvvl          ! variable volume 
    3434   USE trd_oce         ! trends: ocean variables 
    35    USE trdtra          ! trends manager: tracers  
     35   USE trdtra          ! trends manager: tracers 
    3636   USE traqsr          ! penetrative solar radiation (needed for nksr) 
    3737   USE phycst          ! physical constant 
     
    6969      !!                   ***  ROUTINE traatf  *** 
    7070      !! 
    71       !! ** Purpose :   Apply the boundary condition on the after temperature   
     71      !! ** Purpose :   Apply the boundary condition on the after temperature 
    7272      !!             and salinity fields and add the Asselin time filter on now fields. 
    73       !!  
    74       !! ** Method  :   At this stage of the computation, ta and sa are the  
     73      !! 
     74      !! ** Method  :   At this stage of the computation, ta and sa are the 
    7575      !!             after temperature and salinity as the time stepping has 
    7676      !!             been performed in trazdf_imp or trazdf_exp module. 
    7777      !! 
    78       !!              - Apply lateral boundary conditions on (ta,sa)  
    79       !!             at the local domain   boundaries through lbc_lnk call,  
    80       !!             at the one-way open boundaries (ln_bdy=T),  
     78      !!              - Apply lateral boundary conditions on (ta,sa) 
     79      !!             at the local domain   boundaries through lbc_lnk call, 
     80      !!             at the one-way open boundaries (ln_bdy=T), 
    8181      !!             at the AGRIF zoom   boundaries (lk_agrif=T) 
    8282      !! 
     
    8888      INTEGER                                  , INTENT(in   ) :: kt             ! ocean time-step index 
    8989      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Kaa  ! time level indices 
    90       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers  
     90      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers 
    9191      !! 
    9292      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     
    104104 
    105105      ! Update after tracer on domain lateral boundaries 
    106       !  
     106      ! 
    107107#if defined key_agrif 
    108108      CALL Agrif_tra                     ! AGRIF zoom boundaries 
     
    112112      ! 
    113113      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries 
    114   
     114 
    115115      ! set time step size (Euler/Leapfrog) 
    116116      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler) 
     
    119119 
    120120      ! trends computation initialisation 
    121       IF( l_trdtra )   THEN                     
     121      IF( l_trdtra )   THEN 
    122122         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    123123         ztrdt(:,:,jpk) = 0._wp 
    124124         ztrds(:,:,jpk) = 0._wp 
    125          IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend  
     125         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
    126126            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 
    127127            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds ) 
    128128         ENDIF 
    129          ! total trend for the non-time-filtered variables.  
     129         ! total trend for the non-time-filtered variables. 
    130130         zfact = 1.0 / rdt 
    131131         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms 
     
    137137         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds ) 
    138138         IF( ln_linssh ) THEN       ! linear sea surface height only 
    139             ! Store now fields before applying the Asselin filter  
     139            ! Store now fields before applying the Asselin filter 
    140140            ! in order to calculate Asselin filter trend later. 
    141             ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm)  
     141            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 
    142142            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm) 
    143143         ENDIF 
    144144      ENDIF 
    145145 
    146       IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping  
     146      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping 
    147147         ! 
    148148         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl 
     
    156156      ELSE                                            ! Leap-Frog + Asselin filter time stepping 
    157157         ! 
    158          IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface  
     158         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface 
    159159         ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface 
    160160         ENDIF 
     
    164164                  &                    pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1.  ) 
    165165         ! 
    166       ENDIF      
    167       ! 
    168       IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt      
    169          zfact = 1._wp / r2dt              
     166      ENDIF 
     167      ! 
     168      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt 
     169         zfact = 1._wp / r2dt 
    170170         DO jk = 1, jpkm1 
    171171            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact 
     
    191191      !! 
    192192      !! ** Purpose :   fixed volume: apply the Asselin time filter to the "now" field 
    193       !!  
     193      !! 
    194194      !! ** Method  : - Apply a Asselin time filter on now fields. 
    195195      !! 
     
    216216         ! 
    217217         DO_3D_00_00( 1, jpkm1 ) 
    218             ztn = pt(ji,jj,jk,jn,Kmm)                                     
     218            ztn = pt(ji,jj,jk,jn,Kmm) 
    219219            ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers 
    220220            ! 
     
    231231      !!                   ***  ROUTINE tra_atf_vvl  *** 
    232232      !! 
    233       !! ** Purpose :   Time varying volume: apply the Asselin time filter   
    234       !!  
     233      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
     234      !! 
    235235      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
    236236      !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] ) 
     
    262262      ENDIF 
    263263      ! 
    264       IF( cdtype == 'TRA' )  THEN    
     264      IF( cdtype == 'TRA' )  THEN 
    265265         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration 
    266266         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs 
     
    268268      ELSE                          ! passive tracers case 
    269269         ll_traqsr  = .FALSE.          ! NO solar penetration 
    270          ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?   
    271          ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??  
     270         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ? 
     271         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ?? 
    272272      ENDIF 
    273273      ! 
     
    279279      zfact1 = atfp * p2dt 
    280280      zfact2 = zfact1 * r1_rau0 
    281       DO jn = 1, kjpt       
     281      DO jn = 1, kjpt 
    282282         DO_3D_00_00( 1, jpkm1 ) 
    283283            ze3t_b = e3t(ji,jj,jk,Kbb) 
     
    296296            ! 
    297297            ! Add asselin correction on scale factors: 
    298             zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )  
    299             ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) )  
    300             IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * (    rnf_b(ji,jj) -    rnf(ji,jj) )  
     298            zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 
     299            ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) ) 
     300            IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * (    rnf_b(ji,jj) -    rnf(ji,jj) ) 
    301301            IF ( ll_isf ) THEN 
    302302               IF ( ln_isfcav_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) ) 
     
    304304            ENDIF 
    305305            ! 
    306             IF( jk == mikt(ji,jj) ) THEN           ! first level  
     306            IF( jk == mikt(ji,jj) ) THEN           ! first level 
    307307               ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 
    308308            ENDIF 
    309309            ! 
    310310            ! solar penetration (temperature only) 
    311             IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            &  
    312                &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )  
     311            IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
     312               &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
    313313               ! 
    314314            ! 
    315315            IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          & 
    316                &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) &  
     316               &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
    317317               &                              * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) 
    318318 
     
    328328                        &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) 
    329329                  END IF 
    330                   ! level partially include in Losch_2008 ice shelf boundary layer  
     330                  ! level partially include in Losch_2008 ice shelf boundary layer 
    331331                  IF ( jk == misfkb_cav(ji,jj) ) THEN 
    332332                     ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) )  & 
     
    342342                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) 
    343343                  END IF 
    344                   ! level partially include in Losch_2008 ice shelf boundary layer  
     344                  ! level partially include in Losch_2008 ice shelf boundary layer 
    345345                  IF ( jk == misfkb_par(ji,jj) ) THEN 
    346346                     ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  & 
     
    371371            ! 
    372372         END_3D 
    373          !  
     373         ! 
    374374      END DO 
    375375      ! 
    376376      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN 
    377          IF( l_trdtra .AND. cdtype == 'TRA' ) THEN  
     377         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 
    378378            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) ) 
    379379            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) ) 
  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90

    r12481 r12581  
    1 MODULE traatf 
     1MODULE traatfLF 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  traatf  *** 
     3   !!                       ***  MODULE  traatfLF  *** 
    44   !! Ocean active tracers:  Asselin time filtering for temperature and salinity 
    55   !!====================================================================== 
     
    1717   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA 
    1818   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
    19    !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename tranxt.F90 -> traatf.F90. Now only does time filtering. 
     19   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename tranxt.F90 -> traatfLF.F90. Now only does time filtering. 
    2020   !!---------------------------------------------------------------------- 
    2121 
     
    2626   !!---------------------------------------------------------------------- 
    2727   USE oce             ! ocean dynamics and tracers variables 
    28    USE dom_oce         ! ocean space and time domain variables  
     28   USE dom_oce         ! ocean space and time domain variables 
    2929   USE sbc_oce         ! surface boundary condition: ocean 
    3030   USE sbcrnf          ! river runoffs 
     
    3333   USE domvvl          ! variable volume 
    3434   USE trd_oce         ! trends: ocean variables 
    35    USE trdtra          ! trends manager: tracers  
     35   USE trdtra          ! trends manager: tracers 
    3636   USE traqsr          ! penetrative solar radiation (needed for nksr) 
    3737   USE phycst          ! physical constant 
     
    5252   PRIVATE 
    5353 
    54    PUBLIC   tra_atf       ! routine called by step.F90 
    55    PUBLIC   tra_atf_fix   ! to be used in trcnxt 
    56    PUBLIC   tra_atf_vvl   ! to be used in trcnxt 
     54   PUBLIC   tra_atf_lf       ! routine called by step.F90 
     55   PUBLIC   tra_atf_fix_lf   ! to be used in trcnxt !!st WARNING discrepancy here interpol is used 
     56   PUBLIC   tra_atf_qe_lf    ! to be used in trcnxt !!st WARNING discrepancy here interpol is used 
    5757 
    5858   !! * Substitutions 
     
    6565CONTAINS 
    6666 
    67    SUBROUTINE tra_atf( kt, Kbb, Kmm, Kaa, pts ) 
    68       !!---------------------------------------------------------------------- 
    69       !!                   ***  ROUTINE traatf  *** 
    70       !! 
    71       !! ** Purpose :   Apply the boundary condition on the after temperature   
     67   SUBROUTINE tra_atf_lf( kt, Kbb, Kmm, Kaa, pts ) 
     68      !!---------------------------------------------------------------------- 
     69      !!                   ***  ROUTINE traatfLF  *** 
     70      !! 
     71      !! ** Purpose :   Apply the boundary condition on the after temperature 
    7272      !!             and salinity fields and add the Asselin time filter on now fields. 
    73       !!  
    74       !! ** Method  :   At this stage of the computation, ta and sa are the  
     73      !! 
     74      !! ** Method  :   At this stage of the computation, ta and sa are the 
    7575      !!             after temperature and salinity as the time stepping has 
    7676      !!             been performed in trazdf_imp or trazdf_exp module. 
    7777      !! 
    78       !!              - Apply lateral boundary conditions on (ta,sa)  
    79       !!             at the local domain   boundaries through lbc_lnk call,  
    80       !!             at the one-way open boundaries (ln_bdy=T),  
     78      !!              - Apply lateral boundary conditions on (ta,sa) 
     79      !!             at the local domain   boundaries through lbc_lnk call, 
     80      !!             at the one-way open boundaries (ln_bdy=T), 
    8181      !!             at the AGRIF zoom   boundaries (lk_agrif=T) 
    8282      !! 
     
    8888      INTEGER                                  , INTENT(in   ) :: kt             ! ocean time-step index 
    8989      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Kaa  ! time level indices 
    90       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers  
     90      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers 
    9191      !! 
    9292      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     
    9595      !!---------------------------------------------------------------------- 
    9696      ! 
    97       IF( ln_timing )   CALL timing_start( 'tra_atf') 
     97      IF( ln_timing )   CALL timing_start( 'tra_atf_lf') 
    9898      ! 
    9999      IF( kt == nit000 ) THEN 
    100100         IF(lwp) WRITE(numout,*) 
    101          IF(lwp) WRITE(numout,*) 'tra_atf : apply Asselin time filter to "now" fields' 
     101         IF(lwp) WRITE(numout,*) 'tra_atf_lf : apply Asselin time filter to "now" fields' 
    102102         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    103103      ENDIF 
    104104 
    105105      ! Update after tracer on domain lateral boundaries 
    106       !  
     106      ! 
    107107#if defined key_agrif 
    108108      CALL Agrif_tra                     ! AGRIF zoom boundaries 
    109109#endif 
    110110      !                                              ! local domain boundaries  (T-point, unchanged sign) 
    111       CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 
     111      CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 
    112112      ! 
    113113      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries 
    114   
     114 
    115115      ! set time step size (Euler/Leapfrog) 
    116116      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler) 
     
    119119 
    120120      ! trends computation initialisation 
    121       IF( l_trdtra )   THEN                     
     121      IF( l_trdtra )   THEN 
    122122         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    123123         ztrdt(:,:,jpk) = 0._wp 
    124124         ztrds(:,:,jpk) = 0._wp 
    125          IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend  
     125         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
    126126            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 
    127127            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds ) 
    128128         ENDIF 
    129          ! total trend for the non-time-filtered variables.  
     129         ! total trend for the non-time-filtered variables. 
    130130         zfact = 1.0 / rdt 
    131131         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms 
     
    137137         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds ) 
    138138         IF( ln_linssh ) THEN       ! linear sea surface height only 
    139             ! Store now fields before applying the Asselin filter  
     139            ! Store now fields before applying the Asselin filter 
    140140            ! in order to calculate Asselin filter trend later. 
    141             ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm)  
     141            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 
    142142            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm) 
    143143         ENDIF 
    144144      ENDIF 
    145145 
    146       IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping  
     146      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping 
    147147         ! 
    148148         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl 
     
    156156      ELSE                                            ! Leap-Frog + Asselin filter time stepping 
    157157         ! 
    158          IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface  
    159          ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface 
    160          ENDIF 
    161          ! 
    162          CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., & 
     158         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface 
     159         ELSE                   ;   CALL tra_atf_qe_lf( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface 
     160         ENDIF 
     161         ! 
     162         CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., & 
    163163                  &                    pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., & 
    164164                  &                    pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1.  ) 
    165165         ! 
    166       ENDIF      
    167       ! 
    168       IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt      
    169          zfact = 1._wp / r2dt              
     166      ENDIF 
     167      ! 
     168      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt 
     169         zfact = 1._wp / r2dt 
    170170         DO jk = 1, jpkm1 
    171171            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact 
     
    181181         &                                  tab3d_2=pts(:,:,:,jp_sal,Kmm), clinfo2=       ' Sn: ', mask2=tmask ) 
    182182      ! 
    183       IF( ln_timing )   CALL timing_stop('tra_atf') 
    184       ! 
    185    END SUBROUTINE tra_atf 
    186  
    187  
    188    SUBROUTINE tra_atf_fix( kt, Kbb, Kmm, Kaa, kit000, cdtype, pt, kjpt ) 
     183      IF( ln_timing )   CALL timing_stop('tra_atf_lf') 
     184      ! 
     185   END SUBROUTINE tra_atf_lf 
     186 
     187 
     188   SUBROUTINE tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, kit000, cdtype, pt, kjpt ) 
    189189      !!---------------------------------------------------------------------- 
    190190      !!                   ***  ROUTINE tra_atf_fix  *** 
    191191      !! 
    192192      !! ** Purpose :   fixed volume: apply the Asselin time filter to the "now" field 
    193       !!  
     193      !! 
    194194      !! ** Method  : - Apply a Asselin time filter on now fields. 
    195195      !! 
     
    209209      IF( kt == kit000 )  THEN 
    210210         IF(lwp) WRITE(numout,*) 
    211          IF(lwp) WRITE(numout,*) 'tra_atf_fix : time filtering', cdtype 
     211         IF(lwp) WRITE(numout,*) 'tra_atf_fix_lf : time filtering', cdtype 
    212212         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    213213      ENDIF 
     
    216216         ! 
    217217         DO_3D_00_00( 1, jpkm1 ) 
    218             ztn = pt(ji,jj,jk,jn,Kmm)                                     
     218            ztn = pt(ji,jj,jk,jn,Kmm) 
    219219            ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers 
    220220            ! 
     
    224224      END DO 
    225225      ! 
    226    END SUBROUTINE tra_atf_fix 
    227  
    228  
    229    SUBROUTINE tra_atf_vvl( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt ) 
     226   END SUBROUTINE tra_atf_fix_lf 
     227 
     228 
     229   SUBROUTINE tra_atf_qe_lf( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt ) 
    230230      !!---------------------------------------------------------------------- 
    231231      !!                   ***  ROUTINE tra_atf_vvl  *** 
    232232      !! 
    233       !! ** Purpose :   Time varying volume: apply the Asselin time filter   
    234       !!  
     233      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
     234      !! 
    235235      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
    236236      !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] ) 
     
    252252      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
    253253      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar 
    254       REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d, zscale  !   -      - 
     254      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d  !   -      - 
    255255      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf 
    256256      !!---------------------------------------------------------------------- 
     
    262262      ENDIF 
    263263      ! 
    264       IF( cdtype == 'TRA' )  THEN    
     264      IF( cdtype == 'TRA' )  THEN 
    265265         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration 
    266266         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs 
     
    268268      ELSE                          ! passive tracers case 
    269269         ll_traqsr  = .FALSE.          ! NO solar penetration 
    270          ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?   
    271          ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??  
     270         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ? 
     271         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ?? 
    272272      ENDIF 
    273273      ! 
     
    279279      zfact1 = atfp * p2dt 
    280280      zfact2 = zfact1 * r1_rau0 
    281       DO jn = 1, kjpt       
     281      DO jn = 1, kjpt 
    282282         DO_3D_00_00( 1, jpkm1 ) 
    283283            ze3t_b = e3t(ji,jj,jk,Kbb) 
     
    292292            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b 
    293293            ! 
    294             ze3t_f = ze3t_n + atfp * ze3t_d 
    295294            ztc_f  = ztc_n  + atfp * ztc_d 
    296295            ! 
    297             ! Add asselin correction on scale factors: 
    298             zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )  
    299             ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) )  
    300             IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * (    rnf_b(ji,jj) -    rnf(ji,jj) )  
    301             IF ( ll_isf ) THEN 
    302                IF ( ln_isfcav_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) ) 
    303                IF ( ln_isfpar_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) ) 
    304             ENDIF 
    305             ! 
    306             IF( jk == mikt(ji,jj) ) THEN           ! first level  
     296            ! Asselin correction on scale factors is done via ssh in r3t_f 
     297            ze3t_f = e3t_0(ji,jj,jk) * ( 1._wp + r3t_f(ji,jj) * tmask(ji,jj,jk) ) 
     298 
     299            ! 
     300            IF( jk == mikt(ji,jj) ) THEN           ! first level 
    307301               ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 
    308302            ENDIF 
    309303            ! 
    310304            ! solar penetration (temperature only) 
    311             IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            &  
    312                &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )  
     305            IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
     306               &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
    313307               ! 
    314308            ! 
    315309            IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          & 
    316                &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) &  
     310               &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
    317311               &                              * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) 
    318312 
     
    328322                        &                     * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) 
    329323                  END IF 
    330                   ! level partially include in Losch_2008 ice shelf boundary layer  
     324                  ! level partially include in Losch_2008 ice shelf boundary layer 
    331325                  IF ( jk == misfkb_cav(ji,jj) ) THEN 
    332326                     ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) )  & 
     
    342336                            &                 * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) 
    343337                  END IF 
    344                   ! level partially include in Losch_2008 ice shelf boundary layer  
     338                  ! level partially include in Losch_2008 ice shelf boundary layer 
    345339                  IF ( jk == misfkb_par(ji,jj) ) THEN 
    346340                     ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  & 
     
    356350                     ztc_f  = ztc_f  + zfact1 * risfcpl_tsc(ji,jj,jk,jn) * r1_e1e2t(ji,jj) 
    357351                     ! Shouldn't volume increment be spread according thanks to zscale  ? 
    358                      ze3t_f = ze3t_f - zfact1 * risfcpl_vol(ji,jj,jk   ) * r1_e1e2t(ji,jj) 
    359352                  END IF 
    360353                  ! 
     
    371364            ! 
    372365         END_3D 
    373          !  
     366         ! 
    374367      END DO 
    375368      ! 
    376369      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN 
    377          IF( l_trdtra .AND. cdtype == 'TRA' ) THEN  
     370         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 
    378371            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) ) 
    379372            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) ) 
     
    387380      ENDIF 
    388381      ! 
    389    END SUBROUTINE tra_atf_vvl 
     382   END SUBROUTINE tra_atf_qe_lf 
    390383 
    391384   !!====================================================================== 
    392 END MODULE traatf 
     385END MODULE traatfLF 
  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90

    r12492 r12581  
    4141   USE iom              ! xIOs server 
    4242   USE domqe 
     43   USE traatfLF          ! time filtering                   (tra_atf_lf routine) 
     44   USE dynatfLF          ! time filtering                   (dyn_atf_lf routine) 
    4345 
    4446   IMPLICIT NONE 
     
    286288!! 
    287289!!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 
    288                          CALL tra_atf       ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays 
    289                          CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors 
    290290                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height 
     291                         CALL dom_qe_r3c    ( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f )            ! "now" ssh/h_0 ratio from filtrered ssh 
     292                         CALL tra_atf_lf    ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays 
     293                         CALL dyn_atf_lf    ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors 
    291294      ! 
    292295      ! Swap time levels 
Note: See TracChangeset for help on using the changeset viewer.