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_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP/trcatf.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496
  • Property svn:keywords set to Id
File size: 17.4 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#if defined key_mpi3
105      CALL lbc_lnk_nc_multi( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )   
106#else
107      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )   
108#endif
109
110      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Kaa )
111
112      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
113         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
114         ztrdt(:,:,:,:)  = 0._wp
115         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend
116            DO jn = 1, jptra
117               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
118            ENDDO
119         ENDIF
120
121         ! total trend for the non-time-filtered variables.
122         zfact = 1.0 / rn_Dt
123         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3ta*Ta)/e3tn; e3tn cancel from ts(Kmm) terms
124         IF( ln_linssh ) THEN       ! linear sea surface height only
125            DO jn = 1, jptra
126               DO jk = 1, jpkm1
127                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) * zfact
128               END DO
129            END DO
130         ELSE
131            DO jn = 1, jptra
132               DO jk = 1, jpkm1
133                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact
134               END DO
135            END DO
136         ENDIF
137         !
138         DO jn = 1, jptra
139            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
140         ENDDO
141         !
142         IF( ln_linssh ) THEN       ! linear sea surface height only
143            ! Store now fields before applying the Asselin filter
144            ! in order to calculate Asselin filter trend later.
145            ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm) 
146         ENDIF
147
148      ENDIF
149      !                                ! Leap-Frog + Asselin filter time stepping
150      IF( l_1st_euler .OR. ln_top_euler ) THEN    ! Euler time-stepping
151         !
152         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
153            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
154            ztrdt(:,:,:,:) = 0._wp           
155            DO jn = 1, jptra
156               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
157            ENDDO
158         END IF
159         !
160      ELSE     
161         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
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            ENDIF
170         ELSE
171                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )       ! offline
172         ENDIF
173         !
174#if defined key_mpi3
175         CALL lbc_lnk_nc_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
176#else
177         CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
178#endif
179      ENDIF
180      !
181      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
182         DO jn = 1, jptra
183            DO jk = 1, jpkm1
184               zfact = 1._wp / rDt_trc 
185               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
186            END DO
187            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
188         END DO
189      END IF
190      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
191      !
192      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
193         WRITE(charout, FMT="('nxt')")
194         CALL prt_ctl_info( charout, cdcomp = 'top' )
195         CALL prt_ctl(tab4d_1=ptr(:,:,:,:,Kmm), mask1=tmask, clinfo=ctrcnm)
196      ENDIF
197      !
198      IF( ln_timing )   CALL timing_stop('trc_atf')
199      !
200   END SUBROUTINE trc_atf
201
202# if ! defined key_qco
203   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
204      !!----------------------------------------------------------------------
205      !!                   ***  ROUTINE tra_atf_off  ***
206      !!
207      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
208      !!
209      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
210      !!
211      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
212      !!              - save in (ta,sa) a thickness weighted average over the three
213      !!             time levels which will be used to compute rdn and thus the semi-
214      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
215      !!              - swap tracer fields to prepare the next time_step.
216      !!                This can be summurized for tempearture as:
217      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
218      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
219      !!             ztm = 0                                                       otherwise
220      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
221      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
222      !!             tn  = ta
223      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
224      !!
225      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
226      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
227      !!----------------------------------------------------------------------
228      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
229      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
230      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
231      !!     
232      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
233      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
234      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
235      !!----------------------------------------------------------------------
236      !
237      IF( kt == nittrc000 )  THEN
238         IF(lwp) WRITE(numout,*)
239         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
240         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
241         IF( .NOT. ln_linssh ) THEN
242            rfact1 = rn_atfp * rn_Dt
243            rfact2 = rfact1 / rho0
244         ENDIF
245        ! 
246      ENDIF
247      !
248      DO jn = 1, jptra     
249         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
250            ze3t_b = e3t(ji,jj,jk,Kbb)
251            ze3t_n = e3t(ji,jj,jk,Kmm)
252            ze3t_a = e3t(ji,jj,jk,Kaa)
253            !                                         ! tracer content at Before, now and after
254            ztc_b  = ptr(ji,jj,jk,jn,Kbb)  * ze3t_b
255            ztc_n  = ptr(ji,jj,jk,jn,Kmm)  * ze3t_n
256            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
257            !
258            ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
259            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
260            !
261            ze3t_f = ze3t_n + rn_atfp * ze3t_d
262            ztc_f  = ztc_n  + rn_atfp * ztc_d
263            !
264            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
265               ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
266               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
267            ENDIF
268
269            ze3t_f = 1.e0 / ze3t_f
270            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
271            !
272         END_3D
273         !
274      END DO
275      !
276   END SUBROUTINE trc_atf_off
277# else
278   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
279      !!----------------------------------------------------------------------
280      !!                   ***  ROUTINE tra_atf_off  ***
281      !!
282      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
283      !!
284      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
285      !!
286      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
287      !!              - save in (ta,sa) a thickness weighted average over the three
288      !!             time levels which will be used to compute rdn and thus the semi-
289      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
290      !!              - swap tracer fields to prepare the next time_step.
291      !!                This can be summurized for tempearture as:
292      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
293      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
294      !!             ztm = 0                                                       otherwise
295      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
296      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
297      !!             tn  = ta
298      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
299      !!
300      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
301      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
302      !!----------------------------------------------------------------------
303      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
304      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
305      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
306      !!     
307      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
308      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
309      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f           !   -      -
310      !!----------------------------------------------------------------------
311      !
312      IF( kt == nittrc000 )  THEN
313         IF(lwp) WRITE(numout,*)
314         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
315         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
316         IF( .NOT. ln_linssh ) THEN
317            rfact1 = rn_atfp * rn_Dt
318            rfact2 = rfact1 / rho0
319         ENDIF
320        ! 
321      ENDIF
322      !
323      DO jn = 1, jptra     
324         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
325            ze3t_b = 1._wp + r3t(ji,jj,Kbb) * tmask(ji,jj,jk)
326            ze3t_n = 1._wp + r3t(ji,jj,Kmm) * tmask(ji,jj,jk)
327            ze3t_a = 1._wp + r3t(ji,jj,Kaa) * tmask(ji,jj,jk)
328            !                                         ! tracer content at Before, now and after
329            ztc_b  = ptr(ji,jj,jk,jn,Kbb) * ze3t_b
330            ztc_n  = ptr(ji,jj,jk,jn,Kmm) * ze3t_n
331            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
332            !
333            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
334            !
335            ze3t_f = 1._wp + r3t_f(ji,jj)*tmask(ji,jj,jk)
336            ztc_f  = ztc_n  + rn_atfp * ztc_d
337            !
338            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
339               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
340            ENDIF
341
342            ze3t_f = 1.e0 / ze3t_f
343            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
344            !
345         END_3D
346         !
347      END DO
348      !
349   END SUBROUTINE trc_atf_off
350# endif
351
352#else
353   !!----------------------------------------------------------------------
354   !!   Default option                                         Empty module
355   !!----------------------------------------------------------------------
356   USE par_oce
357   USE par_trc
358CONTAINS
359   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr ) 
360      INTEGER                                   , INTENT(in)    :: kt
361      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm, Kaa ! time level indices
362      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr           ! passive tracers and RHS of tracer equation
363      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
364   END SUBROUTINE trc_atf
365#endif
366   !!======================================================================
367END MODULE trcatf
Note: See TracBrowser for help on using the repository browser.