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_r12558_HPC-08_epico_Extra_Halo/src/TOP/TRP – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/TOP/TRP/trcatf.F90 @ 13248

Last change on this file since 13248 was 13248, checked in by francesca, 4 years ago

dev_r12558_HPC-08_epico_Extra_Halo: merge with trunk@13237, see #2366
  • Property svn:keywords set to Id
File size: 17.2 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          ! 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            ENDDO
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         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.
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            ENDDO
154         END IF
155         !
156      ELSE     
157         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
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
162            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
163            ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rn_Dt, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
164# endif
165            ENDIF
166         ELSE
167                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )       ! offline
168         ENDIF
169         !
170         CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
171      ENDIF
172      !
173      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
174         DO jn = 1, jptra
175            DO jk = 1, jpkm1
176               zfact = 1._wp / rDt_trc 
177               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
178            END DO
179            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
180         END DO
181      END IF
182      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
183      !
184      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
185         WRITE(charout, FMT="('nxt')")
186         CALL prt_ctl_info( charout, cdcomp = 'top' )
187         CALL prt_ctl(tab4d_1=ptr(:,:,:,:,Kmm), mask1=tmask, clinfo=ctrcnm)
188      ENDIF
189      !
190      IF( ln_timing )   CALL timing_stop('trc_atf')
191      !
192   END SUBROUTINE trc_atf
193
194# if ! defined key_qco
195   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
196      !!----------------------------------------------------------------------
197      !!                   ***  ROUTINE tra_atf_off  ***
198      !!
199      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
200      !!
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
210      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
211      !!             ztm = 0                                                       otherwise
212      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
213      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
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      !!----------------------------------------------------------------------
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
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,*)
231         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
232         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
233         IF( .NOT. ln_linssh ) THEN
234            rfact1 = rn_atfp * rn_Dt
235            rfact2 = rfact1 / rho0
236         ENDIF
237        ! 
238      ENDIF
239      !
240      DO jn = 1, jptra     
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            !
253            ze3t_f = ze3t_n + rn_atfp * ze3t_d
254            ztc_f  = ztc_n  + rn_atfp * ztc_d
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
260
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
265         !
266      END DO
267      !
268   END SUBROUTINE trc_atf_off
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
333
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
344#else
345   !!----------------------------------------------------------------------
346   !!   Default option                                         Empty module
347   !!----------------------------------------------------------------------
348   USE par_oce
349   USE par_trc
350CONTAINS
351   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr ) 
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
355      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
356   END SUBROUTINE trc_atf
357#endif
358   !!======================================================================
359END MODULE trcatf
Note: See TracBrowser for help on using the repository browser.