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_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/TOP/TRP – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/TOP/TRP/trcatf.F90 @ 14219

Last change on this file since 14219 was 14219, checked in by mcastril, 4 years ago

Add Mixed Precision support by Oriol Tintó
  • 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 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#  include "single_precision_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
61   !! $Id$
62   !! Software governed by the CeCILL license (see ./LICENSE)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )
67      !!----------------------------------------------------------------------
68      !!                   ***  ROUTINE trcatf  ***
69      !!
70      !! ** Purpose :   Apply the boundary condition on the after passive tracers fields and
71      !!      apply Asselin time filter to the now passive tracer fields if using leapfrog timestep
72      !!
73      !! ** Method  :   Apply lateral boundary conditions on (uu(Kaa),vv(Kaa)) through
74      !!      call to lbc_lnk routine
75      !!
76      !!   For Arakawa or TVD Scheme :
77      !!      A Asselin time filter applied on now tracers tr(Kmm) to avoid
78      !!      the divergence of two consecutive time-steps and tr arrays
79      !!      to prepare the next time_step:
80      !!         (tr(Kmm)) = (tr(Kmm)) + rn_atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
81      !!
82      !!
83      !! ** Action  : - update tr(Kmm), tr(Kaa)
84      !!----------------------------------------------------------------------
85      INTEGER                                   , INTENT( in )  :: kt             ! ocean time-step index
86      INTEGER                                   , INTENT( in )  :: Kbb, Kmm, Kaa ! time level indices
87      REAL(dp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers
88      !
89      INTEGER  ::   jk, jn   ! dummy loop indices
90      REAL(wp) ::   zfact            ! temporary scalar
91      CHARACTER (len=22) :: charout
92      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace
93      !!----------------------------------------------------------------------
94      !
95      IF( ln_timing )   CALL timing_start('trc_atf')
96      !
97      IF( kt == nittrc000 .AND. lwp ) THEN
98         WRITE(numout,*)
99         WRITE(numout,*) 'trc_atf : Asselin time filtering on passive tracers'
100      ENDIF
101      !
102#if defined key_agrif
103      CALL Agrif_trc                   ! AGRIF zoom boundaries
104#endif
105      ! Update after tracer on domain lateral boundaries
106      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1._wp )   
107
108      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Kaa )
109
110      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
111         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
112         ztrdt(:,:,:,:)  = 0._wp
113         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend
114            DO jn = 1, jptra
115               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
116            ENDDO
117         ENDIF
118
119         ! total trend for the non-time-filtered variables.
120         zfact = 1.0 / rn_Dt
121         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3ta*Ta)/e3tn; e3tn cancel from ts(Kmm) terms
122         IF( ln_linssh ) THEN       ! linear sea surface height only
123            DO jn = 1, jptra
124               DO jk = 1, jpkm1
125                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) * zfact
126               END DO
127            END DO
128         ELSE
129            DO jn = 1, jptra
130               DO jk = 1, jpkm1
131                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact
132               END DO
133            END DO
134         ENDIF
135         !
136         DO jn = 1, jptra
137            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
138         ENDDO
139         !
140         IF( ln_linssh ) THEN       ! linear sea surface height only
141            ! Store now fields before applying the Asselin filter
142            ! in order to calculate Asselin filter trend later.
143            ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm) 
144         ENDIF
145
146      ENDIF
147      !                                ! Leap-Frog + Asselin filter time stepping
148      IF( l_1st_euler .OR. ln_top_euler ) THEN    ! Euler time-stepping
149         !
150         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
151            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
152            ztrdt(:,:,:,:) = 0._wp           
153            DO jn = 1, jptra
154               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
155            ENDDO
156         END IF
157         !
158      ELSE     
159         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
160# if defined key_qco   ||   defined key_linssh
161            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nittrc000,        'TRC', ptr, jptra )                     !     linear ssh
162            ELSE                   ;   CALL tra_atf_qco_lf( kt, Kbb, Kmm, Kaa, nittrc000, rn_Dt, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
163# else
164            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                       !     linear ssh
165            ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nittrc000, CASTWP(rn_Dt), 'TRC', ptr, sbc_trc, sbc_trc_b, jptra )    ! non-linear ssh
166# endif
167            ENDIF
168         ELSE
169                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )       ! offline
170         ENDIF
171         !
172         CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp )
173      ENDIF
174      !
175      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
176         DO jn = 1, jptra
177            DO jk = 1, jpkm1
178               zfact = 1._wp / rDt_trc 
179               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
180            END DO
181            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
182         END DO
183      END IF
184      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
185      !
186      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
187         WRITE(charout, FMT="('nxt')")
188         CALL prt_ctl_info( charout, cdcomp = 'top' )
189         CALL prt_ctl(tab4d_1=CASTWP(ptr(:,:,:,:,Kmm)), mask1=tmask, clinfo=ctrcnm)
190      ENDIF
191      !
192      IF( ln_timing )   CALL timing_stop('trc_atf')
193      !
194   END SUBROUTINE trc_atf
195
196# if defined key_qco   ||   defined key_linssh
197   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
198      !!----------------------------------------------------------------------
199      !!                   ***  ROUTINE tra_atf_off  ***
200      !!
201      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
202      !!
203      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
204      !!
205      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
206      !!              - save in (ta,sa) a thickness weighted average over the three
207      !!             time levels which will be used to compute rdn and thus the semi-
208      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
209      !!              - swap tracer fields to prepare the next time_step.
210      !!                This can be summurized for tempearture as:
211      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
212      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
213      !!             ztm = 0                                                       otherwise
214      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
215      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
216      !!             tn  = ta
217      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
218      !!
219      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
220      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
221      !!----------------------------------------------------------------------
222      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
223      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
224      REAL(dp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
225      !!     
226      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
227      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
228      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f           !   -      -
229      !!----------------------------------------------------------------------
230      !
231      IF( kt == nittrc000 )  THEN
232         IF(lwp) WRITE(numout,*)
233         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
234         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
235         IF( .NOT. ln_linssh ) THEN
236            rfact1 = rn_atfp * rn_Dt
237            rfact2 = rfact1 / rho0
238         ENDIF
239        ! 
240      ENDIF
241      !
242      DO jn = 1, jptra     
243         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
244            ze3t_b = 1._wp + r3t(ji,jj,Kbb) * tmask(ji,jj,jk)
245            ze3t_n = 1._wp + r3t(ji,jj,Kmm) * tmask(ji,jj,jk)
246            ze3t_a = 1._wp + r3t(ji,jj,Kaa) * tmask(ji,jj,jk)
247            !                                         ! tracer content at Before, now and after
248            ztc_b  = ptr(ji,jj,jk,jn,Kbb) * ze3t_b
249            ztc_n  = ptr(ji,jj,jk,jn,Kmm) * ze3t_n
250            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
251            !
252            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
253            !
254            ze3t_f = 1._wp + r3t_f(ji,jj)*tmask(ji,jj,jk)
255            ztc_f  = ztc_n  + rn_atfp * ztc_d
256            !
257            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
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(dp), 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, ze3t_d   !   -      -
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( 1, 1, 1, 1, 1, jpkm1 )
317            ze3t_b = e3t(ji,jj,jk,Kbb)
318            ze3t_n = e3t(ji,jj,jk,Kmm)
319            ze3t_a = e3t(ji,jj,jk,Kaa)
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            ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
326            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
327            !
328            ze3t_f = ze3t_n + rn_atfp * ze3t_d
329            ztc_f  = ztc_n  + rn_atfp * ztc_d
330            !
331            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
332               ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
333               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
334            ENDIF
335
336            ze3t_f = 1.e0 / ze3t_f
337            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
338            !
339         END_3D
340         !
341      END DO
342      !
343   END SUBROUTINE trc_atf_off
344# endif
345
346#else
347   !!----------------------------------------------------------------------
348   !!   Default option                                         Empty module
349   !!----------------------------------------------------------------------
350   USE par_oce
351   USE par_trc
352CONTAINS
353   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr ) 
354      INTEGER                                   , INTENT(in)    :: kt
355      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm, Kaa ! time level indices
356      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr           ! passive tracers and RHS of tracer equation
357      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
358   END SUBROUTINE trc_atf
359#endif
360   !!======================================================================
361END MODULE trcatf
Note: See TracBrowser for help on using the repository browser.