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_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP – NEMO

source: NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcatf.F90 @ 13373

Last change on this file since 13373 was 13373, checked in by cetlod, 4 years ago

TOP-05_Ethe_Agrif : 1st step of changes to successfully compile, see ticket #2508

  • 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 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
35   USE traatfqco
36# else
37   USE traatf
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                   ! AGRIF zoom boundaries
103#endif
104      ! Update after tracer on domain lateral boundaries
105      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )   
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
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_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), '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
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, ze3t_d   !   -      -
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 = e3t(ji,jj,jk,Kbb)
244            ze3t_n = e3t(ji,jj,jk,Kmm)
245            ze3t_a = e3t(ji,jj,jk,Kaa)
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            ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
252            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
253            !
254            ze3t_f = ze3t_n + rn_atfp * ze3t_d
255            ztc_f  = ztc_n  + rn_atfp * ztc_d
256            !
257            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
258               ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
259               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
260            ENDIF
261
262            ze3t_f = 1.e0 / ze3t_f
263            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
264            !
265         END_3D
266         !
267      END DO
268      !
269   END SUBROUTINE trc_atf_off
270# else
271   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
272      !!----------------------------------------------------------------------
273      !!                   ***  ROUTINE tra_atf_off  ***
274      !!
275      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
276      !!
277      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
278      !!
279      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
280      !!              - save in (ta,sa) a thickness weighted average over the three
281      !!             time levels which will be used to compute rdn and thus the semi-
282      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
283      !!              - swap tracer fields to prepare the next time_step.
284      !!                This can be summurized for tempearture as:
285      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
286      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
287      !!             ztm = 0                                                       otherwise
288      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
289      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
290      !!             tn  = ta
291      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
292      !!
293      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
294      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
295      !!----------------------------------------------------------------------
296      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
297      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
298      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
299      !!     
300      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
301      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
302      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f           !   -      -
303      !!----------------------------------------------------------------------
304      !
305      IF( kt == nittrc000 )  THEN
306         IF(lwp) WRITE(numout,*)
307         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
308         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
309         IF( .NOT. ln_linssh ) THEN
310            rfact1 = rn_atfp * rn_Dt
311            rfact2 = rfact1 / rho0
312         ENDIF
313       
314      ENDIF
315      !
316      DO jn = 1, jptra     
317         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
318            ze3t_b = 1._wp + r3t(ji,jj,Kbb) * tmask(ji,jj,jk)
319            ze3t_n = 1._wp + r3t(ji,jj,Kmm) * tmask(ji,jj,jk)
320            ze3t_a = 1._wp + r3t(ji,jj,Kaa) * tmask(ji,jj,jk)
321            !                                         ! tracer content at Before, now and after
322            ztc_b  = ptr(ji,jj,jk,jn,Kbb) * ze3t_b
323            ztc_n  = ptr(ji,jj,jk,jn,Kmm) * ze3t_n
324            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
325            !
326            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
327            !
328            ze3t_f = 1._wp + r3t_f(ji,jj)*tmask(ji,jj,jk)
329            ztc_f  = ztc_n  + rn_atfp * ztc_d
330            !
331            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
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.