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.
tranxt.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 8698

Last change on this file since 8698 was 8698, checked in by davestorkey, 6 years ago

Bug fixes for tracer trends diagnostics in trunk. See ticket #1877.
These fixes correspond to the changes to the 3.6 branch at revisions 8102 and 8305.
  • Property svn:keywords set to Id
File size: 18.6 KB
Line 
1MODULE tranxt
2   !!======================================================================
3   !!                       ***  MODULE  tranxt  ***
4   !! Ocean active tracers:  time stepping on temperature and salinity
5   !!======================================================================
6   !! History :  OPA  !  1991-11  (G. Madec)  Original code
7   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
8   !!            8.0  !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
9   !!             -   !  1996-04  (A. Weaver)  Euler forward step
10   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
11   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
12   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
13   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget
14   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
15   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
16   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
17   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
18   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   tra_nxt       : time stepping on tracers
23   !!   tra_nxt_fix   : time stepping on tracers : fixed    volume case
24   !!   tra_nxt_vvl   : time stepping on tracers : variable volume case
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers variables
27   USE dom_oce         ! ocean space and time domain variables
28   USE sbc_oce         ! surface boundary condition: ocean
29   USE sbcrnf          ! river runoffs
30   USE sbcisf          ! ice shelf melting
31   USE zdf_oce         ! ocean vertical mixing
32   USE domvvl          ! variable volume
33   USE trd_oce         ! trends: ocean variables
34   USE trdtra          ! trends manager: tracers
35   USE traqsr          ! penetrative solar radiation (needed for nksr)
36   USE phycst          ! physical constant
37   USE ldftra          ! lateral physics on tracers
38   USE ldfslp
39   USE bdy_oce   , ONLY: ln_bdy
40   USE bdytra          ! open boundary condition (bdy_tra routine)
41   !
42   USE in_out_manager  ! I/O manager
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE prtctl          ! Print control
45   USE wrk_nemo        ! Memory allocation
46   USE timing          ! Timing
47#if defined key_agrif
48   USE agrif_opa_interp
49#endif
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC   tra_nxt       ! routine called by step.F90
55   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
56   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
57
58   !! * Substitutions
59#  include "vectopt_loop_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE tra_nxt( kt )
68      !!----------------------------------------------------------------------
69      !!                   ***  ROUTINE tranxt  ***
70      !!
71      !! ** Purpose :   Apply the boundary condition on the after temperature 
72      !!             and salinity fields, achieved the time stepping by adding
73      !!             the Asselin filter on now fields and swapping the fields.
74      !!
75      !! ** Method  :   At this stage of the computation, ta and sa are the
76      !!             after temperature and salinity as the time stepping has
77      !!             been performed in trazdf_imp or trazdf_exp module.
78      !!
79      !!              - Apply lateral boundary conditions on (ta,sa)
80      !!             at the local domain   boundaries through lbc_lnk call,
81      !!             at the one-way open boundaries (ln_bdy=T),
82      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
83      !!
84      !!              - Update lateral boundary conditions on AGRIF children
85      !!             domains (lk_agrif=T)
86      !!
87      !! ** Action  : - tsb & tsn ready for the next time step
88      !!----------------------------------------------------------------------
89      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
90      !!
91      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
92      REAL(wp) ::   zfact            ! local scalars
93      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
94      !!----------------------------------------------------------------------
95      !
96      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt')
97      !
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
101         IF(lwp) WRITE(numout,*) '~~~~~~~'
102      ENDIF
103
104      ! Update after tracer on domain lateral boundaries
105      !
106#if defined key_agrif
107      CALL Agrif_tra                     ! AGRIF zoom boundaries
108#endif
109      !
110      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign)
111      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
112      !
113      IF( ln_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
114 
115      ! set time step size (Euler/Leapfrog)
116      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =     rdt      ! at nit000             (Euler)
117      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
118      ENDIF
119
120      ! trends computation initialisation
121      IF( l_trdtra )   THEN                   
122         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
123         ztrdt(:,:,jpk) = 0._wp
124         ztrds(:,:,jpk) = 0._wp
125         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
126            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
127            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
128         ENDIF
129         ! total trend for the non-time-filtered variables.
130         zfact = 1.0 / rdt
131         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms
132         DO jk = 1, jpkm1
133            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jp_tem)) * zfact
134            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jp_sal)) * zfact
135         END DO
136         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
137         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
138         IF( ln_linssh ) THEN 
139            ! Store now fields before applying the Asselin filter
140            ! in order to calculate Asselin filter trend later.
141            ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
142            ztrds(:,:,:) = tsn(:,:,:,jp_sal)
143         ENDIF
144      ENDIF
145
146      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
147         DO jn = 1, jpts
148            DO jk = 1, jpkm1
149               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
150            END DO
151         END DO
152         IF (l_trdtra .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            ztrds(:,:,:) = 0._wp
156            CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
157            CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
158         END IF
159         !
160      ELSE                                            ! Leap-Frog + Asselin filter time stepping
161         !
162         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface
163         ELSE                   ;   CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa,   &
164           &                                                                sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
165         ENDIF
166         !
167         DO jn = 1, jpts
168            CALL lbc_lnk( tsb(:,:,:,jn), 'T', 1._wp ) 
169            CALL lbc_lnk( tsn(:,:,:,jn), 'T', 1._wp )
170            CALL lbc_lnk( tsa(:,:,:,jn), 'T', 1._wp )
171         END DO
172      ENDIF     
173      !
174      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
175         zfact = 1._wp / r2dt             
176         DO jk = 1, jpkm1
177            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
178            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
179         END DO
180         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
181         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
182      END IF
183      IF( l_trdtra ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
184      !
185      !                        ! control print
186      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
187         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
188      !
189      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
190      !
191   END SUBROUTINE tra_nxt
192
193
194   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
195      !!----------------------------------------------------------------------
196      !!                   ***  ROUTINE tra_nxt_fix  ***
197      !!
198      !! ** Purpose :   fixed volume: apply the Asselin time filter and
199      !!                swap the tracer fields.
200      !!
201      !! ** Method  : - Apply a Asselin time filter on now fields.
202      !!              - swap tracer fields to prepare the next time_step.
203      !!
204      !! ** Action  : - tsb & tsn ready for the next time step
205      !!----------------------------------------------------------------------
206      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
207      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
208      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
209      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
210      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
211      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
212      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
213      !
214      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
215      REAL(wp) ::   ztn, ztd         ! local scalars
216      !!----------------------------------------------------------------------
217      !
218      IF( kt == kit000 )  THEN
219         IF(lwp) WRITE(numout,*)
220         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
221         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
222      ENDIF
223      !
224      DO jn = 1, kjpt
225         !
226         DO jk = 1, jpkm1
227            DO jj = 2, jpjm1
228               DO ji = fs_2, fs_jpim1
229                  ztn = ptn(ji,jj,jk,jn)                                   
230                  ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn)  ! time laplacian on tracers
231                  !
232                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                      ! ptb <-- filtered ptn
233                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                      ! ptn <-- pta
234               END DO
235           END DO
236         END DO
237         !
238      END DO
239      !
240   END SUBROUTINE tra_nxt_fix
241
242
243   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
244      !!----------------------------------------------------------------------
245      !!                   ***  ROUTINE tra_nxt_vvl  ***
246      !!
247      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
248      !!                and swap the tracer fields.
249      !!
250      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
251      !!              - swap tracer fields to prepare the next time_step.
252      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
253      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
254      !!             tn  = ta
255      !!
256      !! ** Action  : - tsb & tsn ready for the next time step
257      !!----------------------------------------------------------------------
258      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
259      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
260      REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step
261      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
262      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
263      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
264      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
265      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
266      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc   ! surface tracer content
267      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
268      !
269      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
270      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
271      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
272      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
273      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrd_atf
274      !!----------------------------------------------------------------------
275      !
276      IF( kt == kit000 )  THEN
277         IF(lwp) WRITE(numout,*)
278         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
279         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
280      ENDIF
281      !
282      IF( cdtype == 'TRA' )  THEN   
283         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
284         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
285         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
286      ELSE                          ! passive tracers case
287         ll_traqsr  = .FALSE.          ! NO solar penetration
288         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ? 
289         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??
290      ENDIF
291      !
292      IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) )   THEN
293         CALL wrk_alloc( jpi, jpj, jpk, kjpt, ztrd_atf )
294         ztrd_atf(:,:,:,:) = 0.0_wp
295      ENDIF
296      zfact = 1._wp / r2dt
297      DO jn = 1, kjpt     
298         DO jk = 1, jpkm1
299            zfact1 = atfp * p2dt
300            zfact2 = zfact1 * r1_rau0
301            DO jj = 2, jpjm1
302               DO ji = fs_2, fs_jpim1
303                  ze3t_b = e3t_b(ji,jj,jk)
304                  ze3t_n = e3t_n(ji,jj,jk)
305                  ze3t_a = e3t_a(ji,jj,jk)
306                  !                                         ! tracer content at Before, now and after
307                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
308                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
309                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
310                  !
311                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
312                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
313                  !
314                  ze3t_f = ze3t_n + atfp * ze3t_d
315                  ztc_f  = ztc_n  + atfp * ztc_d
316                  !
317                  IF( jk == mikt(ji,jj) ) THEN           ! first level
318                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
319                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
320                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
321                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
322                  ENDIF
323                  !
324                  ! solar penetration (temperature only)
325                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
326                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
327                     !
328                  ! river runoff
329                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
330                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
331                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
332                     !
333                  ! ice shelf
334                  IF( ll_isf ) THEN
335                     ! level fully include in the Losch_2008 ice shelf boundary layer
336                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
337                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
338                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
339                     ! level partially include in Losch_2008 ice shelf boundary layer
340                     IF ( jk == misfkb(ji,jj) )                                                   &
341                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
342                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
343                  END IF
344                  !
345                  ze3t_f = 1.e0 / ze3t_f
346                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
347                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
348                  !
349                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
350                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
351                  ENDIF
352                  !
353               END DO
354            END DO
355         END DO
356         !
357      END DO
358      !
359      IF( l_trdtra .and. cdtype == 'TRA' ) THEN
360         CALL trd_tra( kt, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
361         CALL trd_tra( kt, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
362         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
363      ENDIF
364      IF( l_trdtrc .and. cdtype == 'TRC' ) THEN
365         DO jn = 1, kjpt
366            CALL trd_tra( kt, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
367         END DO
368         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
369      ENDIF
370      !
371   END SUBROUTINE tra_nxt_vvl
372
373   !!======================================================================
374END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.