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

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 7555

Last change on this file since 7555 was 7555, checked in by davestorkey, 7 years ago

Fix two bugs in tracer trends diagnostics. Ticket #1823.
  • Property svn:keywords set to Id
File size: 19.8 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/freezing
31   USE zdf_oce         ! ocean vertical mixing
32   USE domvvl          ! variable volume
33   USE dynspg_oce      ! surface     pressure gradient variables
34   USE dynhpg          ! hydrostatic pressure gradient
35   USE trd_oce         ! trends: ocean variables
36   USE trdtra          ! trends manager: tracers
37   USE traqsr          ! penetrative solar radiation (needed for nksr)
38   USE phycst          ! physical constant
39   USE ldftra_oce      ! lateral physics on tracers
40   USE bdy_oce         ! BDY open boundary condition variables
41   USE bdytra          ! open boundary condition (bdy_tra routine)
42   !
43   USE in_out_manager  ! I/O manager
44   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
45   USE prtctl          ! Print control
46   USE wrk_nemo        ! Memory allocation
47   USE timing          ! Timing
48#if defined key_agrif
49   USE agrif_opa_interp
50#endif
51
52   IMPLICIT NONE
53   PRIVATE
54
55   PUBLIC   tra_nxt       ! routine called by step.F90
56   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
57   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
58
59   REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg
60
61   !! * Substitutions
62#  include "domzgr_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
65   !! $Id$
66   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE tra_nxt( kt )
71      !!----------------------------------------------------------------------
72      !!                   ***  ROUTINE tranxt  ***
73      !!
74      !! ** Purpose :   Apply the boundary condition on the after temperature 
75      !!             and salinity fields, achieved the time stepping by adding
76      !!             the Asselin filter on now fields and swapping the fields.
77      !!
78      !! ** Method  :   At this stage of the computation, ta and sa are the
79      !!             after temperature and salinity as the time stepping has
80      !!             been performed in trazdf_imp or trazdf_exp module.
81      !!
82      !!              - Apply lateral boundary conditions on (ta,sa)
83      !!             at the local domain   boundaries through lbc_lnk call,
84      !!             at the one-way open boundaries (lk_bdy=T),
85      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
86      !!
87      !!              - Update lateral boundary conditions on AGRIF children
88      !!             domains (lk_agrif=T)
89      !!
90      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
91      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
92      !!----------------------------------------------------------------------
93      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
94      !!
95      INTEGER  ::   jk, jn    ! dummy loop indices
96      REAL(wp) ::   zfact     ! local scalars
97      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
98      !!----------------------------------------------------------------------
99      !
100      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt')
101      !
102      IF( kt == nit000 ) THEN
103         IF(lwp) WRITE(numout,*)
104         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
105         IF(lwp) WRITE(numout,*) '~~~~~~~'
106         !
107         rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp)      ! Brown & Campana parameter for semi-implicit hpg
108      ENDIF
109
110      ! Update after tracer on domain lateral boundaries
111      !
112#if defined key_agrif
113      CALL Agrif_tra                     ! AGRIF zoom boundaries
114#endif
115      !
116      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign)
117      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
118      !
119#if defined key_bdy 
120      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
121#endif
122 
123      ! set time step size (Euler/Leapfrog)
124      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
125      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
126      ENDIF
127
128      ! trends computation initialisation
129      IF( l_trdtra )   THEN                   
130         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
131         ztrdt(:,:,jk) = 0._wp
132         ztrds(:,:,jk) = 0._wp
133         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
134            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
135            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
136         ENDIF
137         ! total trend for the non-time-filtered variables.
138         DO jk = 1, jpkm1
139            zfact = 1.0 / rdttra(jk)
140            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsn(:,:,jk,jp_tem) ) * zfact 
141            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsn(:,:,jk,jp_sal) ) * zfact 
142         END DO
143         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
144         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
145         ! Store now fields before applying the Asselin filter
146         ! in order to calculate Asselin filter trend later.
147         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
148         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
149      ENDIF
150
151      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
152         DO jn = 1, jpts
153            DO jk = 1, jpkm1
154               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
155            END DO
156         END DO
157      ELSE                                            ! Leap-Frog + Asselin filter time stepping
158         !
159         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   &
160           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)
161         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
162         ENDIF
163      ENDIF     
164      !
165     ! trends computation
166      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
167         DO jk = 1, jpkm1
168            zfact = 1._wp / r2dtra(jk)             
169            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
170            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
171         END DO
172         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
173         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
174         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
175      END IF
176      !
177      !                        ! control print
178      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
179         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
180      !
181      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
182      !
183   END SUBROUTINE tra_nxt
184
185
186   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
187      !!----------------------------------------------------------------------
188      !!                   ***  ROUTINE tra_nxt_fix  ***
189      !!
190      !! ** Purpose :   fixed volume: apply the Asselin time filter and
191      !!                swap the tracer fields.
192      !!
193      !! ** Method  : - Apply a Asselin time filter on now fields.
194      !!              - save in (ta,sa) an average over the three time levels
195      !!             which will be used to compute rdn and thus the semi-implicit
196      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
197      !!              - swap tracer fields to prepare the next time_step.
198      !!                This can be summurized for tempearture as:
199      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
200      !!             ztm = 0                                   otherwise
201      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
202      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
203      !!             tn  = ta 
204      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
205      !!
206      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
207      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
208      !!----------------------------------------------------------------------
209      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
210      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
211      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
212      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
213      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
214      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
215      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
216      !
217      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
218      LOGICAL  ::   ll_tra_hpg       ! local logical
219      REAL(wp) ::   ztn, ztd         ! local scalars
220      !!----------------------------------------------------------------------
221
222      IF( kt == kit000 )  THEN
223         IF(lwp) WRITE(numout,*)
224         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
225         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
226      ENDIF
227      !
228      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
229      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
230      ENDIF
231      !
232      DO jn = 1, kjpt
233         !
234         DO jk = 1, jpkm1
235            DO jj = 1, jpj
236               DO ji = 1, jpi
237                  ztn = ptn(ji,jj,jk,jn)                                   
238                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
239                  !
240                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
241                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
242                  !
243                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
244               END DO
245           END DO
246         END DO
247         !
248      END DO
249      !
250   END SUBROUTINE tra_nxt_fix
251
252
253   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
254      !!----------------------------------------------------------------------
255      !!                   ***  ROUTINE tra_nxt_vvl  ***
256      !!
257      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
258      !!                and swap the tracer fields.
259      !!
260      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
261      !!              - save in (ta,sa) a thickness weighted average over the three
262      !!             time levels which will be used to compute rdn and thus the semi-
263      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
264      !!              - swap tracer fields to prepare the next time_step.
265      !!                This can be summurized for tempearture as:
266      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
267      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
268      !!             ztm = 0                                                       otherwise
269      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
270      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
271      !!             tn  = ta
272      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
273      !!
274      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
275      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
276      !!----------------------------------------------------------------------
277      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
278      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
279      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
280      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
281      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
282      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
283      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
284      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
285      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
286      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
287
288      !!     
289      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf, ll_isf   ! local logical
290      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
291      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
292      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
293      !!----------------------------------------------------------------------
294      !
295      IF( kt == kit000 )  THEN
296         IF(lwp) WRITE(numout,*)
297         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
298         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
299      ENDIF
300      !
301      IF( cdtype == 'TRA' )  THEN   
302         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
303         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
304         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
305         IF (nn_isf .GE. 1) THEN
306            ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing
307         ELSE
308            ll_isf = .FALSE.
309         END IF
310      ELSE                         
311         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
312         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
313         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
314         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing
315      ENDIF
316      !
317      DO jn = 1, kjpt     
318         DO jk = 1, jpkm1
319            zfact1 = atfp * p2dt(jk)
320            zfact2 = zfact1 / rau0
321            DO jj = 1, jpj
322               DO ji = 1, jpi
323                  ze3t_b = fse3t_b(ji,jj,jk)
324                  ze3t_n = fse3t_n(ji,jj,jk)
325                  ze3t_a = fse3t_a(ji,jj,jk)
326                  !                                         ! tracer content at Before, now and after
327                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
328                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
329                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
330                  !
331                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
332                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
333                  !
334                  ze3t_f = ze3t_n + atfp * ze3t_d
335                  ztc_f  = ztc_n  + atfp * ztc_d
336                  !
337                  IF( jk == mikt(ji,jj) ) THEN           ! first level
338                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
339                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
340                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
341                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
342                  ENDIF
343
344                  ! solar penetration (temperature only)
345                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
346                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
347
348                  ! river runoff
349                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
350                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
351                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
352
353                  ! ice shelf
354                  IF( ll_isf ) THEN
355                     ! level fully include in the Losch_2008 ice shelf boundary layer
356                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
357                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
358                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
359                     ! level partially include in Losch_2008 ice shelf boundary layer
360                     IF ( jk == misfkb(ji,jj) )                                                   &
361                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
362                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
363                  END IF
364
365                  ze3t_f = 1.e0 / ze3t_f
366                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
367                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
368                  !
369                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
370                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
371                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
372                  ENDIF
373               END DO
374            END DO
375         END DO
376         !
377      END DO
378      !
379   END SUBROUTINE tra_nxt_vvl
380
381   !!======================================================================
382END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.