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.
trcatf.F90 in NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/TOP/TRP – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/TOP/TRP/trcatf.F90 @ 13151

Last change on this file since 13151 was 13151, checked in by gm, 4 years ago

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