New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11057 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcatf.F90 – NEMO

Ignore:
Timestamp:
2019-05-24T17:36:26+02:00 (5 years 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).
File:
1 moved

Legend:

Unmodified
Added
Removed
  • 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 
Note: See TracChangeset for help on using the changeset viewer.