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 branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 8865

Last change on this file since 8865 was 8865, checked in by deazer, 6 years ago

Moving Changes from CS15mini config into NEMO main src
Updating TEST configs to run wit this version of the code
all sette tests pass at this revision other than AGRIF
Includes updates to dynnxt and tranxt required for 3D rives in WAD case to be conservative.

Next commit will update naming conventions and tidy the code.

  • Property svn:keywords set to Id
File size: 19.2 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                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
320                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
321                  ENDIF
322                  ! Rivers can be not just at the surface must go down to nk_rnd(ji,jj)
323                  IF( ln_rnf_depth ) THEN
324                     IF( mikt(ji,jj) <=jk .and. jk <= nk_rnf(ji,jj)  ) THEN
325                      ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   )  )*(e3t_n(ji,jj,jk)/h_rnf(ji,jj) ) ! as we have sigma can do that here change later
326                     ENDIF
327                  ELSE
328                      IF( jk == mikt(ji,jj) ) THEN           ! first level
329                         ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   ) ) 
330                      ENDIF
331                  ENDIF
332
333                  !
334                  ! solar penetration (temperature only)
335                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
336                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
337                     !
338                  ! river runoff
339                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
340                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
341                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
342                     !
343                  ! ice shelf
344                  IF( ll_isf ) THEN
345                     ! level fully include in the Losch_2008 ice shelf boundary layer
346                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
347                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
348                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
349                     ! level partially include in Losch_2008 ice shelf boundary layer
350                     IF ( jk == misfkb(ji,jj) )                                                   &
351                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
352                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
353                  END IF
354                  !
355                  ze3t_f = 1.e0 / ze3t_f
356                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
357                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
358                  !
359                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
360                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
361                  ENDIF
362                  !
363               END DO
364            END DO
365         END DO
366         !
367      END DO
368      !
369      IF( l_trdtra .and. cdtype == 'TRA' ) THEN
370         CALL trd_tra( kt, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
371         CALL trd_tra( kt, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
372         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
373      ENDIF
374      IF( l_trdtrc .and. cdtype == 'TRC' ) THEN
375         DO jn = 1, kjpt
376            CALL trd_tra( kt, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
377         END DO
378         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
379      ENDIF
380      !
381   END SUBROUTINE tra_nxt_vvl
382
383   !!======================================================================
384END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.