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

source: branches/2016/dev_r6522_SIMPLIF_3/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 6862

Last change on this file since 6862 was 6862, checked in by lovato, 8 years ago

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