Changeset 11057


Ignore:
Timestamp:
2019-05-24T17:36:26+02:00 (16 months ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Rewrite time filtering of tracers and SSH. This version now bit compares with the base code for ORCA1 starting from a restart.
To do:

  1. Check that it bit compares with an Euler timestep.
  2. Full SETTE tests.
  3. Check that it compiles with C1D.
  4. Check that TOP works (see 1).
Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap
Files:
6 edited
2 moved

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/cfgs/SHARED/field_def_nemo-oce.xml

    r11050 r11057  
    6262        <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                          /> 
    6363        <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                          /> 
    64          <field id="mldzint_1"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
    65          <field id="mldzint_2"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
    66          <field id="mldzint_3"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
    67          <field id="mldzint_4"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
    68          <field id="mldzint_5"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
    69          <field id="mldhtc_1"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    70          <field id="mldhtc_2"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    71          <field id="mldhtc_3"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    72          <field id="mldhtc_4"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    73          <field id="mldhtc_5"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    7464        <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    7565        <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="1e-3*kg/m2" /> 
     
    8777        <!-- variables available with MLE --> 
    8878        <field id="Lf_NHpf"      long_name="MLE: Lf = N H / f"   unit="m" /> 
     79 
     80        <!-- next variables available with ln_zad_Aimp=.true. --> 
     81        <field id="Courant"      long_name="Courant number"                                                                 unit="#"   grid_ref="grid_T_3D" /> 
     82        <field id="wimp"         long_name="Implicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" /> 
     83        <field id="wexp"         long_name="Explicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" /> 
     84        <field id="wi_cff"       long_name="Fraction of implicit vertical velocity"                                         unit="#"   grid_ref="grid_T_3D" /> 
    8985 
    9086        <!-- next variables available with key_diahth --> 
     
    355351        <field id="uoce"         long_name="ocean current along i-axis"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D" /> 
    356352        <field id="uoce_e3u"     long_name="ocean current along i-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_U_3D"  > uoce * e3u </field> 
    357         <field id="uoce2_e3u"    long_name="ocean current along i-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_U_3D"  > uoce * uoce * e3u </field> 
    358353        <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             /> 
    359354        <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"                             /> 
     
    410405        <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" /> 
    411406        <field id="voce_e3v"     long_name="ocean current along j-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_V_3D"  > voce * e3v </field> 
    412         <field id="voce2_e3v"    long_name="ocean current along j-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_V_3D"  > voce * voce * e3v </field> 
    413407        <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             /> 
    414408        <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"                             /> 
     
    981975 
    982976    <field_group id="25h_grid_W" grid_ref="grid_W_3D" operation="instant"> 
    983       <field id="vomecrtz25h"         name="k current 25h mean"                 unit="m/s"      /> 
     977      <field id="vovecrtz25h"         name="k current 25h mean"                 unit="m/s"      /> 
    984978      <field id="avt25h"              name="vertical diffusivity25h mean"       unit="m2/s" /> 
    985979      <field id="avm25h"              name="vertical viscosity 25h mean"        unit="m2/s" /> 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/C1D/step_c1d.F90

    r11001 r11057  
    107107      IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs  )  ! horizontal & vertical advection 
    108108      IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes 
    109                         CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   ) ! vertical mixing 
     109                        CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   )         ! vertical mixing 
    110110                        CALL eos( ts(:,:,:,:,Nnn), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl 
    111       IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   ) ! applied non penetrative convective adjustment on (t,s) 
    112                         CALL tra_nxt( kstp, Nbb, Nnn, Nrhs,     Naa   ) ! tracer fields at next time step 
     111      IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   )         ! applied non penetrative convective adjustment on (t,s) 
     112                        CALL tra_atf( kstp, Nbb, Nnn, Nrhs,     Naa, ts   )     ! time filtering of "now" tracer fields 
    113113 
    114114 
     
    124124      IF( ln_zdfosm  )  CALL dyn_osm    ( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes 
    125125                        CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion 
    126                         CALL dyn_nxt    ( kstp, Nbb, Nnn, Naa )                 ! lateral velocity at next time step 
    127       IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp, Nbb, Nnn, Naa )                 ! swap of sea surface height 
    128  
    129       IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )              ! swap of vertical scale factors 
     126                        CALL dyn_atf    ( kstp, Nbb, Nnn, Naa , uu, vv, e3t, e3u, e3v )  ! time filtering of "now" fields 
     127      IF(.NOT.ln_linssh)CALL ssh_atf    ( kstp, Nbb, Nnn, Naa , ssh )                    ! time filtering of "now" sea surface height 
     128      ! 
     129      ! Swap time levels 
     130      Nrhs = Nbb 
     131      Nbb = Nnn 
     132      Nnn = Naa 
     133      Naa = Nrhs 
     134      ! 
     135      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )                       ! swap of vertical scale factors 
    130136 
    131137      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/sshwzv.F90

    r11053 r11057  
    1313   !!---------------------------------------------------------------------- 
    1414   !!   ssh_nxt       : after ssh 
    15    !!   ssh_swp       : filter ans swap the ssh arrays 
     15   !!   ssh_atf       : filter ans swap the ssh arrays 
    1616   !!   wzv           : compute now vertical velocity 
    1717   !!---------------------------------------------------------------------- 
     
    4343   PUBLIC   wzv        ! called by step.F90 
    4444   PUBLIC   wAimp      ! called by step.F90 
    45    PUBLIC   ssh_swp    ! called by step.F90 
     45   PUBLIC   ssh_atf    ! called by step.F90 
    4646 
    4747   !! * Substitutions 
     
    218218 
    219219 
    220    SUBROUTINE ssh_swp( kt, Kbb, Kmm, Kaa ) 
    221       !!---------------------------------------------------------------------- 
    222       !!                    ***  ROUTINE ssh_nxt  *** 
     220   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh ) 
     221      !!---------------------------------------------------------------------- 
     222      !!                    ***  ROUTINE ssh_atf  *** 
     223      !! 
     224      !!      !!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!! 
    223225      !! 
    224226      !! ** Purpose :   achieve the sea surface  height time stepping by  
     
    237239      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    238240      !!---------------------------------------------------------------------- 
    239       INTEGER, INTENT(in) ::   kt              ! ocean time-step index 
    240       INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time-step index 
     241      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index 
     242      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
     243      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field 
    241244      ! 
    242245      REAL(wp) ::   zcoef   ! local scalar 
    243246      !!---------------------------------------------------------------------- 
    244247      ! 
    245       IF( ln_timing )   CALL timing_start('ssh_swp') 
     248      IF( ln_timing )   CALL timing_start('ssh_atf') 
    246249      ! 
    247250      IF( kt == nit000 ) THEN 
    248251         IF(lwp) WRITE(numout,*) 
    249          IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
     252         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height' 
    250253         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    251254      ENDIF 
    252255      !              !==  Euler time-stepping: no filter, just swap  ==! 
    253       IF ( neuler == 0 .AND. kt == nit000 ) THEN 
    254          ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now    <-- after  (before already = now) 
    255          ! 
    256       ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==! 
    257          !                                                  ! before <-- now filtered 
    258          ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) ) 
    259          IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed 
     256      IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN   ! Only do time filtering for leapfrog timesteps 
     257         !                                                  ! filtered "now" field 
     258         ssh(:,:,Kmm) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) ) 
     259         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed 
    260260            zcoef = atfp * rdt * r1_rau0 
    261             ssh(:,:,Kbb) = ssh(:,:,Kbb) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
     261            ssh(:,:,Kmm) = ssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
    262262               &                             -    rnf_b(:,:) + rnf   (:,:)   & 
    263263               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:) 
    264264         ENDIF 
    265          ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now <-- after 
    266       ENDIF 
    267       ! 
    268       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Kbb), clinfo1=' ssh(:,:,Kbb)  - : ', mask1=tmask ) 
    269       ! 
    270       IF( ln_timing )   CALL timing_stop('ssh_swp') 
    271       ! 
    272    END SUBROUTINE ssh_swp 
     265      ENDIF 
     266      ! 
     267      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Kmm), clinfo1=' ssh(:,:,Kmm)  - : ', mask1=tmask ) 
     268      ! 
     269      IF( ln_timing )   CALL timing_stop('ssh_atf') 
     270      ! 
     271   END SUBROUTINE ssh_atf 
    273272 
    274273   SUBROUTINE wAimp( kt, Kmm ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90

    r11056 r11057  
    1 MODULE tranxt 
     1MODULE traatf 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  tranxt  *** 
     3   !!                       ***  MODULE  traatf  *** 
    44   !! Ocean active tracers:  time stepping on temperature and salinity 
    55   !!====================================================================== 
     
    2020 
    2121   !!---------------------------------------------------------------------- 
    22    !!   tra_nxt       : time stepping on tracers 
    23    !!   tra_nxt_fix   : time stepping on tracers : fixed    volume case 
    24    !!   tra_nxt_vvl   : time stepping on tracers : variable volume case 
     22   !!   tra_atf       : time stepping on tracers 
     23   !!   tra_atf_fix   : time stepping on tracers : fixed    volume case 
     24   !!   tra_atf_vvl   : time stepping on tracers : variable volume case 
    2525   !!---------------------------------------------------------------------- 
    2626   USE oce             ! ocean dynamics and tracers variables 
     
    5151   PRIVATE 
    5252 
    53    PUBLIC   tra_nxt       ! routine called by step.F90 
    54    PUBLIC   tra_nxt_fix   ! to be used in trcnxt 
    55    PUBLIC   tra_nxt_vvl   ! to be used in trcnxt 
     53   PUBLIC   tra_atf       ! routine called by step.F90 
     54   PUBLIC   tra_atf_fix   ! to be used in trcnxt 
     55   PUBLIC   tra_atf_vvl   ! to be used in trcnxt 
    5656 
    5757   !! * Substitutions 
     
    6464CONTAINS 
    6565 
    66    SUBROUTINE tra_nxt( kt, Kbb, Kmm, Krhs, Kaa ) 
    67       !!---------------------------------------------------------------------- 
    68       !!                   ***  ROUTINE tranxt  *** 
     66   SUBROUTINE tra_atf( kt, Kbb, Kmm, Kaa, pts ) 
     67      !!---------------------------------------------------------------------- 
     68      !!                   ***  ROUTINE traatf  *** 
     69      !! 
     70      !!       !!!!!!!!!!!!!!!!!  REWRITE HEADER COMMENTS !!!!!!!!!!!!!!! 
    6971      !! 
    7072      !! ** Purpose :   Apply the boundary condition on the after temperature   
     
    8688      !! ** Action  : - ts(Kbb) & ts(Kmm) ready for the next time step 
    8789      !!---------------------------------------------------------------------- 
    88       INTEGER, INTENT(in) ::   kt                   ! ocean time-step index 
    89       INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs, Kaa  ! time level indices 
     90      INTEGER                                  , INTENT(in   ) :: kt             ! ocean time-step index 
     91      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Kaa  ! time level indices 
     92      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers  
    9093      !! 
    9194      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     
    9497      !!---------------------------------------------------------------------- 
    9598      ! 
    96       IF( ln_timing )   CALL timing_start( 'tra_nxt') 
     99      IF( ln_timing )   CALL timing_start( 'tra_atf') 
    97100      ! 
    98101      IF( kt == nit000 ) THEN 
    99102         IF(lwp) WRITE(numout,*) 
    100          IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 
     103         IF(lwp) WRITE(numout,*) 'tra_atf : achieve the time stepping by Asselin filter and array swap' 
    101104         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    102105      ENDIF 
     
    108111#endif 
    109112      !                                              ! local domain boundaries  (T-point, unchanged sign) 
    110       CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. ) 
    111       ! 
    112       !! IMMERSE development : Following the general pattern for the new code we want to pass in the  
    113       !!                       velocities to bdy_dyn as arguments so here we use "ts" even  
    114       !!                       though we haven't converted the tracer names in the rest of tranxt.F90 
    115       !!                       because it will be completely rewritten. DS. 
    116       IF( ln_bdy )   CALL bdy_tra( kt, Kbb, ts, Kaa )  ! BDY open boundaries 
     113      CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 
     114      ! 
     115      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries 
    117116  
    118117      ! set time step size (Euler/Leapfrog) 
     
    127126         ztrds(:,:,jpk) = 0._wp 
    128127         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend  
    129             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 
    130             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdfp, ztrds ) 
     128            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 
     129            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds ) 
    131130         ENDIF 
    132131         ! total trend for the non-time-filtered variables.  
    133132         zfact = 1.0 / rdt 
    134          ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms 
     133         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms 
    135134         DO jk = 1, jpkm1 
    136             ztrdt(:,:,jk) = ( ts(:,:,jk,jp_tem,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - ts(:,:,jk,jp_tem,Kmm)) * zfact 
    137             ztrds(:,:,jk) = ( ts(:,:,jk,jp_sal,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - ts(:,:,jk,jp_sal,Kmm)) * zfact 
    138          END DO 
    139          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_tot, ztrdt ) 
    140          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_tot, ztrds ) 
     135            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_tem,Kmm)) * zfact 
     136            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_sal,Kmm)) * zfact 
     137         END DO 
     138         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_tot, ztrdt ) 
     139         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds ) 
    141140         IF( ln_linssh ) THEN       ! linear sea surface height only 
    142141            ! Store now fields before applying the Asselin filter  
    143142            ! in order to calculate Asselin filter trend later. 
    144             ztrdt(:,:,:) = ts(:,:,:,jp_tem,Kmm)  
    145             ztrds(:,:,:) = ts(:,:,:,jp_sal,Kmm) 
     143            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm)  
     144            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm) 
    146145         ENDIF 
    147146      ENDIF 
     
    150149         DO jn = 1, jpts 
    151150            DO jk = 1, jpkm1 
    152                ts(:,:,jk,jn,Kmm) = ts(:,:,jk,jn,Krhs)     
     151               pts(:,:,jk,jn,Kmm) = pts(:,:,jk,jn,Kaa)     
    153152            END DO 
    154153         END DO 
    155154         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl 
    156             !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step 
     155            !                                        ! Asselin filter is output by tra_atf_vvl that is not called on this time step 
    157156            ztrdt(:,:,:) = 0._wp 
    158157            ztrds(:,:,:) = 0._wp 
    159             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_atf, ztrdt ) 
    160             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_atf, ztrds ) 
     158            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt ) 
     159            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds ) 
    161160         END IF 
    162161         ! 
    163162      ELSE                                            ! Leap-Frog + Asselin filter time stepping 
    164163         ! 
    165          IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt,      Kmm,       nit000,      'TRA',                           & 
    166            &                                                   ts(:,:,:,:,Kbb), ts(:,:,:,:,Kmm), ts(:,:,:,:,Krhs), jpts )  ! linear free surface  
    167          ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Krhs, nit000, rdt, 'TRA',                           & 
    168            &                                                   ts(:,:,:,:,Kbb), ts(:,:,:,:,Kmm), ts(:,:,:,:,Krhs),      & 
    169            &                                                   sbc_tsc        , sbc_tsc_b                        , jpts )  ! non-linear free surface 
    170          ENDIF 
    171          ! 
    172          CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Kbb) , 'T', 1., ts(:,:,:,jp_sal,Kbb) , 'T', 1., & 
    173                   &                    ts(:,:,:,jp_tem,Kmm) , 'T', 1., ts(:,:,:,jp_sal,Kmm) , 'T', 1., & 
    174                   &                    ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1.  ) 
     164         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface  
     165         ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface 
     166         ENDIF 
     167         ! 
     168         CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., & 
     169                  &                    pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., & 
     170                  &                    pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1.  ) 
    175171         ! 
    176172      ENDIF      
     
    179175         zfact = 1._wp / r2dt              
    180176         DO jk = 1, jpkm1 
    181             ztrdt(:,:,jk) = ( ts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact 
    182             ztrds(:,:,jk) = ( ts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact 
    183          END DO 
    184          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_atf, ztrdt ) 
    185          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_atf, ztrds ) 
     177            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact 
     178            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact 
     179         END DO 
     180         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt ) 
     181         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds ) 
    186182      END IF 
    187183      IF( l_trdtra )   DEALLOCATE( ztrdt , ztrds ) 
    188184      ! 
    189185      !                        ! control print 
    190       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' nxt  - Tn: ', mask1=tmask,   & 
    191          &                       tab3d_2=ts(:,:,:,jp_sal,Kmm), clinfo2=       ' Sn: ', mask2=tmask ) 
    192       ! 
    193       IF( ln_timing )   CALL timing_stop('tra_nxt') 
    194       ! 
    195    END SUBROUTINE tra_nxt 
    196  
    197  
    198    SUBROUTINE tra_nxt_fix( kt, Kmm, kit000, cdtype, ptb, ptn, pta, kjpt ) 
    199       !!---------------------------------------------------------------------- 
    200       !!                   ***  ROUTINE tra_nxt_fix  *** 
     186      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kmm), clinfo1=' nxt  - Tn: ', mask1=tmask,   & 
     187         &                       tab3d_2=pts(:,:,:,jp_sal,Kmm), clinfo2=       ' Sn: ', mask2=tmask ) 
     188      ! 
     189      IF( ln_timing )   CALL timing_stop('tra_atf') 
     190      ! 
     191   END SUBROUTINE tra_atf 
     192 
     193 
     194   SUBROUTINE tra_atf_fix( kt, kit000, Kbb, Kmm, Kaa, cdtype, pt, kjpt ) 
     195      !!---------------------------------------------------------------------- 
     196      !!                   ***  ROUTINE tra_atf_fix  *** 
    201197      !! 
    202198      !! ** Purpose :   fixed volume: apply the Asselin time filter and  
     
    208204      !! ** Action  : - ptb & ptn ready for the next time step 
    209205      !!---------------------------------------------------------------------- 
    210       INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index 
    211       INTEGER                              , INTENT(in   ) ::  Kmm       ! time level index 
    212       INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index 
    213       CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator) 
    214       INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers 
    215       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields 
    216       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields 
    217       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend 
     206      INTEGER                                  , INTENT(in   ) ::  kt            ! ocean time-step index 
     207      INTEGER                                  , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices 
     208      INTEGER                                  , INTENT(in   ) ::  kit000        ! first time step index 
     209      CHARACTER(len=3)                         , INTENT(in   ) ::  cdtype        ! =TRA or TRC (tracer indicator) 
     210      INTEGER                                  , INTENT(in   ) ::  kjpt          ! number of tracers 
     211      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::  pt            ! tracer fields 
    218212      ! 
    219213      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     
    223217      IF( kt == kit000 )  THEN 
    224218         IF(lwp) WRITE(numout,*) 
    225          IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype 
     219         IF(lwp) WRITE(numout,*) 'tra_atf_fix : time stepping', cdtype 
    226220         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    227221      ENDIF 
     
    232226            DO jj = 2, jpjm1 
    233227               DO ji = fs_2, fs_jpim1 
    234                   ztn = ptn(ji,jj,jk,jn)                                     
    235                   ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn)  ! time laplacian on tracers 
    236                   ! 
    237                   ptb(ji,jj,jk,jn) = ztn + atfp * ztd                      ! ptb <-- filtered ptn  
    238                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                      ! ptn <-- pta 
     228                  ztn = pt(ji,jj,jk,jn,Kmm)                                     
     229                  ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers 
     230                  ! 
     231                  pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd                      ! ptb <-- filtered ptn  
    239232               END DO 
    240233           END DO 
     
    243236      END DO 
    244237      ! 
    245    END SUBROUTINE tra_nxt_fix 
    246  
    247  
    248    SUBROUTINE tra_nxt_vvl( kt, Kbb, Kmm, Krhs, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt ) 
    249       !!---------------------------------------------------------------------- 
    250       !!                   ***  ROUTINE tra_nxt_vvl  *** 
     238   END SUBROUTINE tra_atf_fix 
     239 
     240 
     241   SUBROUTINE tra_atf_vvl( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt ) 
     242      !!---------------------------------------------------------------------- 
     243      !!                   ***  ROUTINE tra_atf_vvl  *** 
    251244      !! 
    252245      !! ** Purpose :   Time varying volume: apply the Asselin time filter   
     
    256249      !!              - swap tracer fields to prepare the next time_step. 
    257250      !!             tb  = ( e3t(Kmm)*tn + atfp*[ e3t(Kbb)*tb - 2 e3t(Kmm)*tn + e3t_a*ta ] ) 
    258       !!                  /( e3t(Kmm)    + atfp*[ e3t(Kbb)    - 2 e3t(Kmm)    + e3t(Krhs)    ] ) 
     251      !!                  /( e3t(Kmm)    + atfp*[ e3t(Kbb)    - 2 e3t(Kmm)    + e3t(Kaa)    ] ) 
    259252      !!             tn  = ta  
    260253      !! 
    261254      !! ** Action  : - ptb & ptn ready for the next time step 
    262255      !!---------------------------------------------------------------------- 
    263       INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index 
    264       INTEGER                              , INTENT(in   ) ::  Kbb, Kmm, Krhs ! time level indices 
    265       INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index 
    266       REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step 
    267       CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator) 
    268       INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers 
    269       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields 
    270       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields 
    271       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend 
    272       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc   ! surface tracer content 
    273       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc_b ! before surface tracer content 
     256      INTEGER                                  , INTENT(in   ) ::  kt        ! ocean time-step index 
     257      INTEGER                                  , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices 
     258      INTEGER                                  , INTENT(in   ) ::  kit000    ! first time step index 
     259      REAL(wp)                                 , INTENT(in   ) ::  p2dt      ! time-step 
     260      CHARACTER(len=3)                         , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator) 
     261      INTEGER                                  , INTENT(in   ) ::  kjpt      ! number of tracers 
     262      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::  pt        ! tracer fields 
     263      REAL(wp), DIMENSION(jpi,jpj    ,kjpt)    , INTENT(in   ) ::  psbc_tc   ! surface tracer content 
     264      REAL(wp), DIMENSION(jpi,jpj    ,kjpt)    , INTENT(in   ) ::  psbc_tc_b ! before surface tracer content 
    274265      ! 
    275266      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical 
     
    282273      IF( kt == kit000 )  THEN 
    283274         IF(lwp) WRITE(numout,*) 
    284          IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype 
     275         IF(lwp) WRITE(numout,*) 'tra_atf_vvl : time stepping', cdtype 
    285276         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    286277      ENDIF 
     
    309300                  ze3t_b = e3t(ji,jj,jk,Kbb) 
    310301                  ze3t_n = e3t(ji,jj,jk,Kmm) 
    311                   ze3t_a = e3t(ji,jj,jk,Krhs) 
     302                  ze3t_a = e3t(ji,jj,jk,Kaa) 
    312303                  !                                         ! tracer content at Before, now and after 
    313                   ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
    314                   ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n 
    315                   ztc_a  = pta(ji,jj,jk,jn) * ze3t_a 
     304                  ztc_b  = pt(ji,jj,jk,jn,Kbb) * ze3t_b 
     305                  ztc_n  = pt(ji,jj,jk,jn,Kmm) * ze3t_n 
     306                  ztc_a  = pt(ji,jj,jk,jn,Kaa) * ze3t_a 
    316307                  ! 
    317308                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b 
     
    361352                  ! 
    362353                  ze3t_f = 1.e0 / ze3t_f 
    363                   ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
    364                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
     354                  pt(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f    ! time filtered "now" field 
    365355                  ! 
    366356                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN 
     
    376366      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN 
    377367         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN  
    378             CALL trd_tra( kt, Kmm, Krhs, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) ) 
    379             CALL trd_tra( kt, Kmm, Krhs, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) ) 
     368            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) ) 
     369            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) ) 
    380370         ENDIF 
    381371         IF( l_trdtrc .AND. cdtype == 'TRC' ) THEN 
    382372            DO jn = 1, kjpt 
    383                CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) ) 
     373               CALL trd_tra( kt, Kmm, Kaa, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) ) 
    384374            END DO 
    385375         ENDIF 
     
    387377      ENDIF 
    388378      ! 
    389    END SUBROUTINE tra_nxt_vvl 
     379   END SUBROUTINE tra_atf_vvl 
    390380 
    391381   !!====================================================================== 
    392 END MODULE tranxt 
     382END MODULE traatf 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/step.F90

    r11053 r11057  
    280280!!  
    281281!!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 
    282                          CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! finalize (bcs) tracer fields at next time step and swap 
    283                          CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" arrays 
    284                          CALL ssh_swp       ( kstp, Nbb, Nnn, Naa )  ! swap of sea surface height 
     282                         CALL tra_atf       ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays 
     283                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors 
     284                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height 
    285285      ! 
    286286      ! Swap time levels 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/step_oce.F90

    r11050 r11057  
    3030   USE traldf          ! lateral mixing                   (tra_ldf routine) 
    3131   USE trazdf          ! vertical mixing                  (tra_zdf routine) 
    32    USE tranxt          ! time-stepping                    (tra_nxt routine) 
     32   USE traatf          ! time filtering                   (tra_atf routine) 
    3333   USE tranpc          ! non-penetrative convection       (tra_npc routine) 
    3434 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcatf.F90

    r11056 r11057  
    1 MODULE trcnxt 
     1MODULE trcatf 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  trcnxt  *** 
     3   !!                       ***  MODULE  trcatf  *** 
    44   !! Ocean passive tracers:  time stepping on passives tracers 
    55   !!====================================================================== 
     
    2424   !!   'key_top'                                                TOP models 
    2525   !!---------------------------------------------------------------------- 
    26    !!   trc_nxt     : time stepping on passive tracers 
     26   !!   trc_atf     : time stepping on passive tracers 
    2727   !!---------------------------------------------------------------------- 
    2828   USE oce_trc         ! ocean dynamics and tracers variables 
     
    4343   PRIVATE 
    4444 
    45    PUBLIC   trc_nxt   ! routine called by step.F90 
     45   PUBLIC   trc_atf   ! routine called by step.F90 
    4646 
    4747   REAL(wp) ::   rfact1, rfact2 
     
    5454CONTAINS 
    5555 
    56    SUBROUTINE trc_nxt( kt, Kbb, Kmm, Krhs ) 
    57       !!---------------------------------------------------------------------- 
    58       !!                   ***  ROUTINE trcnxt  *** 
     56   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr ) 
     57      !!---------------------------------------------------------------------- 
     58      !!                   ***  ROUTINE trcatf  *** 
     59      !! 
     60      !!        !!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!! 
    5961      !! 
    6062      !! ** Purpose :   Compute the passive tracers fields at the  
    6163      !!      next time-step from their temporal trends and swap the fields. 
    6264      !!  
    63       !! ** Method  :   Apply lateral boundary conditions on (uu(Krhs),vv(Krhs)) through  
     65      !! ** Method  :   Apply lateral boundary conditions on (uu(Kaa),vv(Kaa)) through  
    6466      !!      call to lbc_lnk routine 
    6567      !!   default: 
    6668      !!      arrays swap 
    67       !!         (tr(Kmm)) = (tr(Krhs)) ; (tr(Krhs)) = (0,0) 
     69      !!         (tr(Kmm)) = (tr(Kaa)) ; (tr(Kaa)) = (0,0) 
    6870      !!         (tr(Kbb)) = (tr(Kmm))  
    6971      !! 
     
    7274      !!      the divergence of two consecutive time-steps and tr arrays 
    7375      !!      to prepare the next time_step: 
    74       !!         (tr(Kbb)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Krhs)) - 2 (tr(Kmm)) ] 
    75       !!         (tr(Kmm)) = (tr(Krhs)) ; (tr(Krhs)) = (0,0) 
     76      !!         (tr(Kbb)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ] 
     77      !!         (tr(Kmm)) = (tr(Kaa)) ; (tr(Kaa)) = (0,0) 
    7678      !! 
    7779      !! 
    7880      !! ** Action  : - update tr(Kbb), tr(Kmm) 
    7981      !!---------------------------------------------------------------------- 
    80       INTEGER, INTENT( in ) ::   kt     ! ocean time-step index 
    81       INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! time level indices 
     82      INTEGER                                   , INTENT( in )  :: kt             ! ocean time-step index 
     83      INTEGER                                   , INTENT( in )  :: Kbb, Kmm, Kaa ! time level indices 
     84      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers 
    8285      ! 
    8386      INTEGER  ::   jk, jn   ! dummy loop indices 
     
    8790      !!---------------------------------------------------------------------- 
    8891      ! 
    89       IF( ln_timing )   CALL timing_start('trc_nxt') 
     92      IF( ln_timing )   CALL timing_start('trc_atf') 
    9093      ! 
    9194      IF( kt == nittrc000 .AND. lwp ) THEN 
    9295         WRITE(numout,*) 
    93          WRITE(numout,*) 'trc_nxt : time stepping on passive tracers' 
     96         WRITE(numout,*) 'trc_atf : Asselin time filtering on passive tracers' 
    9497      ENDIF 
    9598      ! 
     
    98101#endif 
    99102      ! Update after tracer on domain lateral boundaries 
    100       CALL lbc_lnk( 'trcnxt', tr(:,:,:,:,Krhs), 'T', 1. )    
    101  
    102       IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Krhs ) 
     103      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )    
     104 
     105      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Kaa ) 
    103106 
    104107      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application 
     
    107110         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend  
    108111            DO jn = 1, jptra 
    109                CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) ) 
     112               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) ) 
    110113            ENDDO 
    111114         ENDIF 
     
    117120            DO jn = 1, jptra 
    118121               DO jk = 1, jpkm1 
    119                   ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) - tr(:,:,jk,jn,Kmm)) * zfact 
     122                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) * zfact 
    120123               END DO 
    121124            END DO 
     
    123126            DO jn = 1, jptra 
    124127               DO jk = 1, jpkm1 
    125                   ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Krhs) - tr(:,:,jk,jn,Kmm) ) * zfact 
     128                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact 
    126129               END DO 
    127130            END DO 
     
    129132         ! 
    130133         DO jn = 1, jptra 
    131             CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) ) 
     134            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) ) 
    132135         ENDDO 
    133136         ! 
     
    135138            ! Store now fields before applying the Asselin filter  
    136139            ! in order to calculate Asselin filter trend later. 
    137             ztrdt(:,:,:,:) = tr(:,:,:,:,Kmm)  
     140            ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm)  
    138141         ENDIF 
    139142 
    140143      ENDIF 
    141144      !                                ! Leap-Frog + Asselin filter time stepping 
    142       IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping (only swap) 
    143          DO jn = 1, jptra 
    144             DO jk = 1, jpkm1 
    145                tr(:,:,jk,jn,Kmm) = tr(:,:,jk,jn,Krhs) 
    146                tr(:,:,jk,jn,Kbb) = tr(:,:,jk,jn,Kmm)   
    147             END DO 
    148          END DO 
    149          IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl 
    150             !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step 
    151             ztrdt(:,:,:,:) = 0._wp             
    152             DO jn = 1, jptra 
    153                CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) ) 
    154             ENDDO 
    155          END IF 
    156          ! 
    157       ELSE      
     145      IF( .NOT.( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) ) THEN    ! Only time filter if not an Euler timestep 
    158146         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping 
    159             IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt,      Kmm,       nittrc000,         'TRC',                & 
    160               &                                              tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), tr(:,:,:,:,Krhs), jptra )  !     linear ssh 
    161             ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Krhs, nittrc000, rdttrc, 'TRC',                & 
    162               &                                              tr(:,:,:,:,Kbb), tr(:,:,:,:,Kmm), tr(:,:,:,:,Krhs),      & 
    163               &                                              sbc_trc, sbc_trc_b, jptra )                                  ! non-linear ssh 
     147            IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh 
     148            ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rdttrc, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh 
    164149            ENDIF 
    165150         ELSE 
    166                                        CALL trc_nxt_off( kt, Kbb, Kmm, Krhs )       ! offline  
     151                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa )       ! offline  
    167152         ENDIF 
    168153         ! 
    169          CALL lbc_lnk_multi( 'trcnxt', tr(:,:,:,:,Kbb), 'T', 1._wp, tr(:,:,:,:,Kmm), 'T', 1._wp, tr(:,:,:,:,Krhs), 'T', 1._wp ) 
     154         CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kbb), 'T', 1._wp, ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp ) 
    170155      ENDIF 
    171156      ! 
     
    174159            DO jk = 1, jpkm1 
    175160               zfact = 1._wp / r2dttrc   
    176                ztrdt(:,:,jk,jn) = ( tr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact  
    177             END DO 
    178             CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) ) 
     161               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact  
     162            END DO 
     163            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) ) 
    179164         END DO 
    180165      END IF 
     
    184169         WRITE(charout, FMT="('nxt')") 
    185170         CALL prt_ctl_trc_info(charout) 
    186          CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm) 
    187       ENDIF 
    188       ! 
    189       IF( ln_timing )   CALL timing_stop('trc_nxt') 
    190       ! 
    191    END SUBROUTINE trc_nxt 
    192  
    193  
    194    SUBROUTINE trc_nxt_off( kt, Kbb, Kmm, Krhs ) 
    195       !!---------------------------------------------------------------------- 
    196       !!                   ***  ROUTINE tra_nxt_vvl  *** 
     171         CALL prt_ctl_trc(tab4d=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm) 
     172      ENDIF 
     173      ! 
     174      IF( ln_timing )   CALL timing_stop('trc_atf') 
     175      ! 
     176   END SUBROUTINE trc_atf 
     177 
     178 
     179   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr ) 
     180      !!---------------------------------------------------------------------- 
     181      !!                   ***  ROUTINE tra_atf_off  *** 
     182      !! 
     183      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!! 
    197184      !! 
    198185      !! ** Purpose :   Time varying volume: apply the Asselin time filter   
    199       !!                and swap the tracer fields. 
    200186      !!  
    201187      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
     
    206192      !!                This can be summurized for tempearture as: 
    207193      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T 
    208       !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Krhs)    ] )    
     194      !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )    
    209195      !!             ztm = 0                                                       otherwise 
    210196      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
    211       !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Krhs)    ] ) 
     197      !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] ) 
    212198      !!             tn  = ta  
    213199      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call) 
     
    216202      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    217203      !!---------------------------------------------------------------------- 
    218       INTEGER, INTENT(in   )   ::  kt              ! ocean time-step index 
    219       INTEGER, INTENT(in   )   ::  Kbb, Kmm, Krhs  ! time level indices 
     204      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index 
     205      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices 
     206      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers 
    220207      !!      
    221208      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
     
    226213      IF( kt == nittrc000 )  THEN 
    227214         IF(lwp) WRITE(numout,*) 
    228          IF(lwp) WRITE(numout,*) 'trc_nxt_off : time stepping' 
     215         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering' 
    229216         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    230217         IF( .NOT. ln_linssh ) THEN 
     
    241228                  ze3t_b = e3t(ji,jj,jk,Kbb) 
    242229                  ze3t_n = e3t(ji,jj,jk,Kmm) 
    243                   ze3t_a = e3t(ji,jj,jk,Krhs) 
     230                  ze3t_a = e3t(ji,jj,jk,Kaa) 
    244231                  !                                         ! tracer content at Before, now and after 
    245                   ztc_b  = tr(ji,jj,jk,jn,Kbb)  * ze3t_b 
    246                   ztc_n  = tr(ji,jj,jk,jn,Kmm)  * ze3t_n 
    247                   ztc_a  = tr(ji,jj,jk,jn,Krhs) * ze3t_a 
     232                  ztc_b  = ptr(ji,jj,jk,jn,Kbb)  * ze3t_b 
     233                  ztc_n  = ptr(ji,jj,jk,jn,Kmm)  * ze3t_n 
     234                  ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a 
    248235                  ! 
    249236                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b 
     
    259246 
    260247                  ze3t_f = 1.e0 / ze3t_f 
    261                   tr(ji,jj,jk,jn,Kbb) = ztc_f * ze3t_f           ! pt(:,:,:,:,Kbb) <-- pt(:,:,:,:,Kmm) filtered 
    262                   tr(ji,jj,jk,jn,Kmm) = tr(ji,jj,jk,jn,Krhs)     ! pt(:,:,:,:,Kmm) <-- pt(:,:,:,:,Krhs) 
     248                  ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field 
    263249                  ! 
    264250               END DO 
     
    268254      END DO 
    269255      ! 
    270    END SUBROUTINE trc_nxt_off 
     256   END SUBROUTINE trc_atf_off 
    271257 
    272258#else 
     
    275261   !!---------------------------------------------------------------------- 
    276262CONTAINS 
    277    SUBROUTINE trc_nxt( kt, Kbb, Kmm, Krhs )   
    278       INTEGER, INTENT(in) :: kt 
    279       INTEGER, INTENT(in) :: Kbb, Kmm, Krhs  ! time level indices 
    280       WRITE(*,*) 'trc_nxt: You should not have seen this print! error?', kt 
    281    END SUBROUTINE trc_nxt 
     263   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )   
     264      INTEGER                                   , INTENT(in) :: kt 
     265      INTEGER                                   , INTENT(in) :: Kbb, Kmm, Kaa  ! time level indices 
     266      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers 
     267      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt 
     268   END SUBROUTINE trc_atf 
    282269#endif 
    283270   !!====================================================================== 
    284 END MODULE trcnxt 
     271END MODULE trcatf 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trctrp.F90

    r10985 r11057  
    2020   USE trcadv          ! advection                           (trc_adv routine) 
    2121   USE trczdf          ! vertical diffusion                  (trc_zdf routine) 
    22    USE trcnxt          ! time-stepping                       (trc_nxt routine) 
     22   USE trcatf          ! time filtering                      (trc_atf routine) 
    2323   USE trcrad          ! positivity                          (trc_rad routine) 
    2424   USE trcsbc          ! surface boundary condition          (trc_sbc routine) 
     
    7878#endif 
    7979                                CALL trc_zdf    ( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer   ==> after 
    80                                 CALL trc_nxt    ( kt, Kbb, Kmm, Krhs )            ! tracer fields at next time step      
     80                                CALL trc_atf    ( kt, Kbb, Kmm, Kaa , tr )        ! time filtering of "now" tracer fields     
    8181         IF( ln_trcrad )        CALL trc_rad    ( kt, Kbb, Kmm, Krhs, tr       )  ! Correct artificial negative concentrations 
    8282         IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kt, Kbb, Kmm )                  ! internal damping trends on closed seas only 
     
    8787         IF( ln_trcdmp )        CALL trc_dmp( kt, Kbb, Kmm,       tr, Krhs )  ! internal damping trends 
    8888                                CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer ==> after 
    89                                 CALL trc_nxt( kt, Kbb, Kmm, Krhs )            ! tracer fields at next time step      
     89                                CALL trc_atf( kt, Kbb, Kmm, Kaa , tr )        ! time filtering of "now" tracer fields 
    9090          IF( ln_trcrad )       CALL trc_rad( kt, Kbb, Kmm, Krhs, tr       )  ! Correct artificial negative concentrations 
    9191         ! 
Note: See TracChangeset for help on using the changeset viewer.