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/2021/dev_r14318_RK3_stage1/src/TOP/TRP – NEMO

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/TRP/trcatf.F90 @ 14800

Last change on this file since 14800 was 14800, checked in by jchanut, 3 years ago

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