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

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

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 13.0 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   USE traatf
34   USE bdy_oce   , ONLY: ln_bdy
35   USE trcbdy          ! BDY open boundaries
36# if defined key_agrif
37   USE agrif_top_interp
38# endif
39   !
40   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
41   USE prtctl_trc      ! Print control for debbuging
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   trc_atf   ! routine called by step.F90
47
48   REAL(wp) ::   rfact1, rfact2
49
50   !! * Substitutions
51#  include "do_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
54   !! $Id$
55   !! Software governed by the CeCILL license (see ./LICENSE)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )
60      !!----------------------------------------------------------------------
61      !!                   ***  ROUTINE trcatf  ***
62      !!
63      !! ** Purpose :   Apply the boundary condition on the after passive tracers fields and
64      !!      apply Asselin time filter to the now passive tracer fields if using leapfrog timestep
65      !!
66      !! ** Method  :   Apply lateral boundary conditions on (uu(Kaa),vv(Kaa)) through
67      !!      call to lbc_lnk routine
68      !!
69      !!   For Arakawa or TVD Scheme :
70      !!      A Asselin time filter applied on now tracers tr(Kmm) to avoid
71      !!      the divergence of two consecutive time-steps and tr arrays
72      !!      to prepare the next time_step:
73      !!         (tr(Kmm)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
74      !!
75      !!
76      !! ** Action  : - update tr(Kmm), tr(Kaa)
77      !!----------------------------------------------------------------------
78      INTEGER                                   , INTENT( in )  :: kt             ! ocean time-step index
79      INTEGER                                   , INTENT( in )  :: Kbb, Kmm, Kaa ! time level indices
80      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers
81      !
82      INTEGER  ::   jk, jn   ! dummy loop indices
83      REAL(wp) ::   zfact            ! temporary scalar
84      CHARACTER (len=22) :: charout
85      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace
86      !!----------------------------------------------------------------------
87      !
88      IF( ln_timing )   CALL timing_start('trc_atf')
89      !
90      IF( kt == nittrc000 .AND. lwp ) THEN
91         WRITE(numout,*)
92         WRITE(numout,*) 'trc_atf : Asselin time filtering on passive tracers'
93      ENDIF
94      !
95#if defined key_agrif
96      CALL Agrif_trc                   ! AGRIF zoom boundaries
97#endif
98      ! Update after tracer on domain lateral boundaries
99      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )   
100
101      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Kaa )
102
103      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
104         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
105         ztrdt(:,:,:,:)  = 0._wp
106         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend
107            DO jn = 1, jptra
108               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
109            ENDDO
110         ENDIF
111
112         ! total trend for the non-time-filtered variables.
113         zfact = 1.0 / rdttrc
114         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms
115         IF( ln_linssh ) THEN       ! linear sea surface height only
116            DO jn = 1, jptra
117               DO jk = 1, jpkm1
118                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) * zfact
119               END DO
120            END DO
121         ELSE
122            DO jn = 1, jptra
123               DO jk = 1, jpkm1
124                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact
125               END DO
126            END DO
127         ENDIF
128         !
129         DO jn = 1, jptra
130            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
131         ENDDO
132         !
133         IF( ln_linssh ) THEN       ! linear sea surface height only
134            ! Store now fields before applying the Asselin filter
135            ! in order to calculate Asselin filter trend later.
136            ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm) 
137         ENDIF
138
139      ENDIF
140      !                                ! Leap-Frog + Asselin filter time stepping
141      IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping
142         !
143         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
144            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
145            ztrdt(:,:,:,:) = 0._wp           
146            DO jn = 1, jptra
147               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
148            ENDDO
149         END IF
150         !
151      ELSE     
152         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
153            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
154            ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rdttrc, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
155            ENDIF
156         ELSE
157                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )       ! offline
158         ENDIF
159         !
160         CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
161      ENDIF
162      !
163      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
164         DO jn = 1, jptra
165            DO jk = 1, jpkm1
166               zfact = 1._wp / r2dttrc 
167               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
168            END DO
169            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
170         END DO
171      END IF
172      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
173      !
174      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
175         WRITE(charout, FMT="('nxt')")
176         CALL prt_ctl_trc_info(charout)
177         CALL prt_ctl_trc(tab4d=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
178      ENDIF
179      !
180      IF( ln_timing )   CALL timing_stop('trc_atf')
181      !
182   END SUBROUTINE trc_atf
183
184
185   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
186      !!----------------------------------------------------------------------
187      !!                   ***  ROUTINE tra_atf_off  ***
188      !!
189      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
190      !!
191      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
192      !!
193      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
194      !!              - save in (ta,sa) a thickness weighted average over the three
195      !!             time levels which will be used to compute rdn and thus the semi-
196      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
197      !!              - swap tracer fields to prepare the next time_step.
198      !!                This can be summurized for tempearture as:
199      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
200      !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )   
201      !!             ztm = 0                                                       otherwise
202      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
203      !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )
204      !!             tn  = ta
205      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
206      !!
207      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
208      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
209      !!----------------------------------------------------------------------
210      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
211      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
212      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
213      !!     
214      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
215      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
216      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
217      !!----------------------------------------------------------------------
218      !
219      IF( kt == nittrc000 )  THEN
220         IF(lwp) WRITE(numout,*)
221         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
222         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
223         IF( .NOT. ln_linssh ) THEN
224            rfact1 = atfp * rdttrc
225            rfact2 = rfact1 / rau0
226         ENDIF
227        ! 
228      ENDIF
229      !
230      DO jn = 1, jptra     
231         DO_3D_11_11( 1, jpkm1 )
232            ze3t_b = e3t(ji,jj,jk,Kbb)
233            ze3t_n = e3t(ji,jj,jk,Kmm)
234            ze3t_a = e3t(ji,jj,jk,Kaa)
235            !                                         ! tracer content at Before, now and after
236            ztc_b  = ptr(ji,jj,jk,jn,Kbb)  * ze3t_b
237            ztc_n  = ptr(ji,jj,jk,jn,Kmm)  * ze3t_n
238            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
239            !
240            ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
241            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
242            !
243            ze3t_f = ze3t_n + atfp * ze3t_d
244            ztc_f  = ztc_n  + atfp * ztc_d
245            !
246            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
247               ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
248               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
249            ENDIF
250
251            ze3t_f = 1.e0 / ze3t_f
252            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
253            !
254         END_3D
255         !
256      END DO
257      !
258   END SUBROUTINE trc_atf_off
259
260#else
261   !!----------------------------------------------------------------------
262   !!   Default option                                         Empty module
263   !!----------------------------------------------------------------------
264   USE par_oce
265   USE par_trc
266CONTAINS
267   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr ) 
268      INTEGER                                   , INTENT(in)    :: kt
269      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm, Kaa ! time level indices
270      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr           ! passive tracers and RHS of tracer equation
271      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
272   END SUBROUTINE trc_atf
273#endif
274   !!======================================================================
275END MODULE trcatf
Note: See TracBrowser for help on using the repository browser.