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/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 8877

Last change on this file since 8877 was 8877, checked in by frrh, 6 years ago

Clear out SVN keywords and properties.

File size: 16.9 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(:,:,jk) = 0._wp
124         ztrds(:,:,jk) = 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         DO jk = 1, jpkm1
132            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsn(:,:,jk,jp_tem) ) * zfact 
133            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsn(:,:,jk,jp_sal) ) * zfact 
134         END DO
135         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
136         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
137         ! Store now fields before applying the Asselin filter
138         ! in order to calculate Asselin filter trend later.
139         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
140         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
141      ENDIF
142
143      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
144         DO jn = 1, jpts
145            DO jk = 1, jpkm1
146               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
147            END DO
148         END DO
149         !
150      ELSE                                            ! Leap-Frog + Asselin filter time stepping
151         !
152         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface
153         ELSE                   ;   CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa,   &
154           &                                                                sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
155         ENDIF
156         !
157         DO jn = 1, jpts
158            CALL lbc_lnk( tsb(:,:,:,jn), 'T', 1._wp ) 
159            CALL lbc_lnk( tsn(:,:,:,jn), 'T', 1._wp )
160            CALL lbc_lnk( tsa(:,:,:,jn), 'T', 1._wp )
161         END DO
162      ENDIF     
163      !
164      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
165         DO jk = 1, jpkm1
166            zfact = 1._wp / r2dt             
167            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
168            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
169         END DO
170         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
171         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
172         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
173      END IF
174      !
175      !                        ! control print
176      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
177         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
178      !
179      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
180      !
181   END SUBROUTINE tra_nxt
182
183
184   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
185      !!----------------------------------------------------------------------
186      !!                   ***  ROUTINE tra_nxt_fix  ***
187      !!
188      !! ** Purpose :   fixed volume: apply the Asselin time filter and
189      !!                swap the tracer fields.
190      !!
191      !! ** Method  : - Apply a Asselin time filter on now fields.
192      !!              - swap tracer fields to prepare the next time_step.
193      !!
194      !! ** Action  : - tsb & tsn ready for the next time step
195      !!----------------------------------------------------------------------
196      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
197      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
198      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
199      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
200      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
201      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
202      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
203      !
204      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
205      REAL(wp) ::   ztn, ztd         ! local scalars
206      !!----------------------------------------------------------------------
207      !
208      IF( kt == kit000 )  THEN
209         IF(lwp) WRITE(numout,*)
210         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
211         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
212      ENDIF
213      !
214      DO jn = 1, kjpt
215         !
216         DO jk = 1, jpkm1
217            DO jj = 2, jpjm1
218               DO ji = fs_2, fs_jpim1
219                  ztn = ptn(ji,jj,jk,jn)                                   
220                  ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn)  ! time laplacian on tracers
221                  !
222                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                      ! ptb <-- filtered ptn
223                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                      ! ptn <-- pta
224               END DO
225           END DO
226         END DO
227         !
228      END DO
229      !
230   END SUBROUTINE tra_nxt_fix
231
232
233   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
234      !!----------------------------------------------------------------------
235      !!                   ***  ROUTINE tra_nxt_vvl  ***
236      !!
237      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
238      !!                and swap the tracer fields.
239      !!
240      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
241      !!              - swap tracer fields to prepare the next time_step.
242      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
243      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
244      !!             tn  = ta
245      !!
246      !! ** Action  : - tsb & tsn ready for the next time step
247      !!----------------------------------------------------------------------
248      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
249      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
250      REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step
251      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
252      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
253      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
254      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
255      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
256      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc   ! surface tracer content
257      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
258      !
259      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
260      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
261      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
262      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
263      !!----------------------------------------------------------------------
264      !
265      IF( kt == kit000 )  THEN
266         IF(lwp) WRITE(numout,*)
267         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
268         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
269      ENDIF
270      !
271      IF( cdtype == 'TRA' )  THEN   
272         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
273         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
274         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
275      ELSE                          ! passive tracers case
276         ll_traqsr  = .FALSE.          ! NO solar penetration
277         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ? 
278         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??
279      ENDIF
280      !
281      DO jn = 1, kjpt     
282         DO jk = 1, jpkm1
283            zfact1 = atfp * p2dt
284            zfact2 = zfact1 * r1_rau0
285            DO jj = 2, jpjm1
286               DO ji = fs_2, fs_jpim1
287                  ze3t_b = e3t_b(ji,jj,jk)
288                  ze3t_n = e3t_n(ji,jj,jk)
289                  ze3t_a = e3t_a(ji,jj,jk)
290                  !                                         ! tracer content at Before, now and after
291                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
292                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
293                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
294                  !
295                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
296                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
297                  !
298                  ze3t_f = ze3t_n + atfp * ze3t_d
299                  ztc_f  = ztc_n  + atfp * ztc_d
300                  !
301                  IF( jk == mikt(ji,jj) ) THEN           ! first level
302                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
303                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
304                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
305                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
306                  ENDIF
307                  !
308                  ! solar penetration (temperature only)
309                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
310                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
311                     !
312                  ! river runoff
313                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
314                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
315                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
316                     !
317                  ! ice shelf
318                  IF( ll_isf ) THEN
319                     ! level fully include in the Losch_2008 ice shelf boundary layer
320                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
321                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
322                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
323                     ! level partially include in Losch_2008 ice shelf boundary layer
324                     IF ( jk == misfkb(ji,jj) )                                                   &
325                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
326                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
327                  END IF
328                  !
329                  ze3t_f = 1.e0 / ze3t_f
330                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
331                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
332                  !
333               END DO
334            END DO
335         END DO
336         !
337      END DO
338      !
339   END SUBROUTINE tra_nxt_vvl
340
341   !!======================================================================
342END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.