source: NEMO/trunk/src/TOP/TRP/trcatf.F90 @ 13237

Last change on this file since 13237 was 13237, checked in by smasson, 4 months ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

  • Property svn:keywords set to Id
File size: 17.2 KB
RevLine 
[11057]1MODULE trcatf
[941]2   !!======================================================================
[11057]3   !!                       ***  MODULE  trcatf  ***
[941]4   !! Ocean passive tracers:  time stepping on passives tracers
5   !!======================================================================
[2528]6   !! History :  7.0  !  1991-11  (G. Madec)  Original code
7   !!                 !  1993-03  (M. Guyon)  symetrical conditions
8   !!                 !  1995-02  (M. Levy)   passive tracers
9   !!                 !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
10   !!            8.0  !  1996-04  (A. Weaver)  Euler forward step
11   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
12   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
13   !!                 !  2002-08  (G. Madec)  F90: Free form and module
14   !!                 !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
15   !!                 !  2004-03  (C. Ethe) passive tracers
16   !!                 !  2007-02  (C. Deltel) Diagnose ML trends for passive tracers
17   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
18   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
19   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
20   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
[11483]21   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename trcnxt.F90 -> trcatf.F90. Now only does time filtering.
[1175]22   !!----------------------------------------------------------------------
[941]23#if defined key_top
24   !!----------------------------------------------------------------------
25   !!   'key_top'                                                TOP models
26   !!----------------------------------------------------------------------
[11057]27   !!   trc_atf     : time stepping on passive tracers
[941]28   !!----------------------------------------------------------------------
29   USE oce_trc         ! ocean dynamics and tracers variables
[2528]30   USE trc             ! ocean passive tracers variables
[4990]31   USE trd_oce
[2528]32   USE trdtra
[13237]33# if defined key_qco
34   USE traatfqco
35# else
[11099]36   USE traatf
[13237]37# endif
[7646]38   USE bdy_oce   , ONLY: ln_bdy
[6140]39   USE trcbdy          ! BDY open boundaries
[2528]40# if defined key_agrif
[941]41   USE agrif_top_interp
[2528]42# endif
[9019]43   !
44   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
45   USE prtctl_trc      ! Print control for debbuging
[941]46
47   IMPLICIT NONE
48   PRIVATE
49
[11057]50   PUBLIC   trc_atf   ! routine called by step.F90
[2528]51
[7881]52   REAL(wp) ::   rfact1, rfact2
[7872]53
[12340]54   !! * Substitutions
55#  include "do_loop_substitute.h90"
[13237]56#  include "domzgr_substitute.h90"
[941]57   !!----------------------------------------------------------------------
[10067]58   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
[7753]59   !! $Id$
[10068]60   !! Software governed by the CeCILL license (see ./LICENSE)
[941]61   !!----------------------------------------------------------------------
62CONTAINS
63
[11057]64   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )
[941]65      !!----------------------------------------------------------------------
[11057]66      !!                   ***  ROUTINE trcatf  ***
[941]67      !!
[11483]68      !! ** Purpose :   Apply the boundary condition on the after passive tracers fields and
69      !!      apply Asselin time filter to the now passive tracer fields if using leapfrog timestep
[941]70      !!
[11057]71      !! ** Method  :   Apply lateral boundary conditions on (uu(Kaa),vv(Kaa)) through
[941]72      !!      call to lbc_lnk routine
73      !!
74      !!   For Arakawa or TVD Scheme :
[11483]75      !!      A Asselin time filter applied on now tracers tr(Kmm) to avoid
[941]76      !!      the divergence of two consecutive time-steps and tr arrays
77      !!      to prepare the next time_step:
[12489]78      !!         (tr(Kmm)) = (tr(Kmm)) + rn_atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
[941]79      !!
80      !!
[11483]81      !! ** Action  : - update tr(Kmm), tr(Kaa)
[941]82      !!----------------------------------------------------------------------
[11057]83      INTEGER                                   , INTENT( in )  :: kt             ! ocean time-step index
84      INTEGER                                   , INTENT( in )  :: Kbb, Kmm, Kaa ! time level indices
85      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers
[2715]86      !
[7753]87      INTEGER  ::   jk, jn   ! dummy loop indices
[941]88      REAL(wp) ::   zfact            ! temporary scalar
89      CHARACTER (len=22) :: charout
[9019]90      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace
[941]91      !!----------------------------------------------------------------------
[3294]92      !
[11057]93      IF( ln_timing )   CALL timing_start('trc_atf')
[3294]94      !
95      IF( kt == nittrc000 .AND. lwp ) THEN
[941]96         WRITE(numout,*)
[11057]97         WRITE(numout,*) 'trc_atf : Asselin time filtering on passive tracers'
[941]98      ENDIF
[6140]99      !
[5656]100#if defined key_agrif
101      CALL Agrif_trc                   ! AGRIF zoom boundaries
102#endif
[9081]103      ! Update after tracer on domain lateral boundaries
[11057]104      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )   
[941]105
[11057]106      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Kaa )
[941]107
[6140]108      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
[9019]109         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
[10096]110         ztrdt(:,:,:,:)  = 0._wp
111         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend
112            DO jn = 1, jptra
[11057]113               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
[10096]114            ENDDO
115         ENDIF
116
117         ! total trend for the non-time-filtered variables.
[12489]118         zfact = 1.0 / rn_Dt
[13237]119         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3ta*Ta)/e3tn; e3tn cancel from ts(Kmm) terms
[10096]120         IF( ln_linssh ) THEN       ! linear sea surface height only
121            DO jn = 1, jptra
122               DO jk = 1, jpkm1
[11057]123                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) * zfact
[10096]124               END DO
125            END DO
126         ELSE
127            DO jn = 1, jptra
128               DO jk = 1, jpkm1
[11057]129                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact
[10096]130               END DO
131            END DO
132         ENDIF
133         !
134         DO jn = 1, jptra
[11057]135            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
[10096]136         ENDDO
137         !
138         IF( ln_linssh ) THEN       ! linear sea surface height only
139            ! Store now fields before applying the Asselin filter
140            ! in order to calculate Asselin filter trend later.
[11057]141            ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm) 
[10096]142         ENDIF
143
[2528]144      ENDIF
[6140]145      !                                ! Leap-Frog + Asselin filter time stepping
[12489]146      IF( l_1st_euler .OR. ln_top_euler ) THEN    ! Euler time-stepping
[11099]147         !
148         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
149            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
150            ztrdt(:,:,:,:) = 0._wp           
151            DO jn = 1, jptra
152               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
153            ENDDO
154         END IF
155         !
156      ELSE     
[7881]157         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
[13237]158# if defined key_qco
159            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nittrc000,        'TRC', ptr, jptra )                     !     linear ssh
160            ELSE                   ;   CALL tra_atf_qco_lf( kt, Kbb, Kmm, Kaa, nittrc000, rn_Dt, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
161# else
[11099]162            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
[12489]163            ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rn_Dt, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
[13237]164# endif
[7872]165            ENDIF
166         ELSE
[11099]167                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )       ! offline
[2528]168         ENDIF
[6140]169         !
[11427]170         CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
[2528]171      ENDIF
[6140]172      !
[10096]173      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
[2528]174         DO jn = 1, jptra
175            DO jk = 1, jpkm1
[12489]176               zfact = 1._wp / rDt_trc 
[11057]177               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
[1175]178            END DO
[11057]179            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
[2528]180         END DO
181      END IF
[10096]182      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
[2528]183      !
[12236]184      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
[941]185         WRITE(charout, FMT="('nxt')")
186         CALL prt_ctl_trc_info(charout)
[11057]187         CALL prt_ctl_trc(tab4d=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
[941]188      ENDIF
[2528]189      !
[11057]190      IF( ln_timing )   CALL timing_stop('trc_atf')
[3294]191      !
[11057]192   END SUBROUTINE trc_atf
[941]193
[13237]194# if ! defined key_qco
[11057]195   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
[7872]196      !!----------------------------------------------------------------------
[11057]197      !!                   ***  ROUTINE tra_atf_off  ***
[7872]198      !!
[11057]199      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
200      !!
[7872]201      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
202      !!
203      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
204      !!              - save in (ta,sa) a thickness weighted average over the three
205      !!             time levels which will be used to compute rdn and thus the semi-
206      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
207      !!              - swap tracer fields to prepare the next time_step.
208      !!                This can be summurized for tempearture as:
209      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
[13237]210      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
[7872]211      !!             ztm = 0                                                       otherwise
[12489]212      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
[13237]213      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
[7872]214      !!             tn  = ta
215      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
216      !!
217      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
218      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
219      !!----------------------------------------------------------------------
[11057]220      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
221      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
222      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
[7872]223      !!     
224      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
225      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
226      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
227      !!----------------------------------------------------------------------
228      !
229      IF( kt == nittrc000 )  THEN
230         IF(lwp) WRITE(numout,*)
[11057]231         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
[7872]232         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
233         IF( .NOT. ln_linssh ) THEN
[12489]234            rfact1 = rn_atfp * rn_Dt
235            rfact2 = rfact1 / rho0
[7872]236         ENDIF
237       
238      ENDIF
239      !
240      DO jn = 1, jptra     
[12340]241         DO_3D_11_11( 1, jpkm1 )
242            ze3t_b = e3t(ji,jj,jk,Kbb)
243            ze3t_n = e3t(ji,jj,jk,Kmm)
244            ze3t_a = e3t(ji,jj,jk,Kaa)
245            !                                         ! tracer content at Before, now and after
246            ztc_b  = ptr(ji,jj,jk,jn,Kbb)  * ze3t_b
247            ztc_n  = ptr(ji,jj,jk,jn,Kmm)  * ze3t_n
248            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
249            !
250            ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
251            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
252            !
[12489]253            ze3t_f = ze3t_n + rn_atfp * ze3t_d
254            ztc_f  = ztc_n  + rn_atfp * ztc_d
[12340]255            !
256            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
257               ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
258               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
259            ENDIF
[7872]260
[12340]261            ze3t_f = 1.e0 / ze3t_f
262            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
263            !
264         END_3D
[7872]265         !
266      END DO
267      !
[11057]268   END SUBROUTINE trc_atf_off
[13237]269# else
270   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
271      !!----------------------------------------------------------------------
272      !!                   ***  ROUTINE tra_atf_off  ***
273      !!
274      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
275      !!
276      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
277      !!
278      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
279      !!              - save in (ta,sa) a thickness weighted average over the three
280      !!             time levels which will be used to compute rdn and thus the semi-
281      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
282      !!              - swap tracer fields to prepare the next time_step.
283      !!                This can be summurized for tempearture as:
284      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
285      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
286      !!             ztm = 0                                                       otherwise
287      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
288      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
289      !!             tn  = ta
290      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
291      !!
292      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
293      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
294      !!----------------------------------------------------------------------
295      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
296      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
297      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
298      !!     
299      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
300      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
301      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f           !   -      -
302      !!----------------------------------------------------------------------
303      !
304      IF( kt == nittrc000 )  THEN
305         IF(lwp) WRITE(numout,*)
306         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
307         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
308         IF( .NOT. ln_linssh ) THEN
309            rfact1 = rn_atfp * rn_Dt
310            rfact2 = rfact1 / rho0
311         ENDIF
312       
313      ENDIF
314      !
315      DO jn = 1, jptra     
316         DO_3D_11_11( 1, jpkm1 )
317            ze3t_b = 1._wp + r3t(ji,jj,Kbb) * tmask(ji,jj,jk)
318            ze3t_n = 1._wp + r3t(ji,jj,Kmm) * tmask(ji,jj,jk)
319            ze3t_a = 1._wp + r3t(ji,jj,Kaa) * tmask(ji,jj,jk)
320            !                                         ! tracer content at Before, now and after
321            ztc_b  = ptr(ji,jj,jk,jn,Kbb) * ze3t_b
322            ztc_n  = ptr(ji,jj,jk,jn,Kmm) * ze3t_n
323            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
324            !
325            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
326            !
327            ze3t_f = 1._wp + r3t_f(ji,jj)*tmask(ji,jj,jk)
328            ztc_f  = ztc_n  + rn_atfp * ztc_d
329            !
330            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
331               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
332            ENDIF
[7872]333
[13237]334            ze3t_f = 1.e0 / ze3t_f
335            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
336            !
337         END_3D
338         !
339      END DO
340      !
341   END SUBROUTINE trc_atf_off
342# endif
343
[941]344#else
345   !!----------------------------------------------------------------------
346   !!   Default option                                         Empty module
347   !!----------------------------------------------------------------------
[11099]348   USE par_oce
349   USE par_trc
[941]350CONTAINS
[11057]351   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr ) 
[11099]352      INTEGER                                   , INTENT(in)    :: kt
353      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm, Kaa ! time level indices
354      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr           ! passive tracers and RHS of tracer equation
[11057]355      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
356   END SUBROUTINE trc_atf
[941]357#endif
358   !!======================================================================
[11057]359END MODULE trcatf
Note: See TracBrowser for help on using the repository browser.