source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 7061

Last change on this file since 7061 was 7061, checked in by davestorkey, 5 years ago

Enable 3D tracer trends diagnostics ,ie. diagnostics obtained for ln_tra_trd=T, in UKMO/dev_r5518_GO6_package branch.
Remove error check for VVL and implement thickness weighting.
Also implement new diagnostics:
ttrd_totad : the total trend due to advection (correct alternative to ttrd_ad!)
ttrd_iso_x/y/z1 : individual components of trend due to isopycnal diffusion. Note the "z1" diagnostic is only part of the vertical trend. The other part is ttrd_zdf minus ttrd_zdfp.
ttrd_evd : the trend due to the EVD convection term (which is contained in ttrd_zdf and ttrd_zdfp).
ttrd_tot : the total model trend.
u/v/weiv_masstr3d : the components of the mass transport due to the Gent-Mc Williams? scheme.
u/v/weiv_heattr3d : the components of the heat transport due to the Gent-Mc Williams? scheme.
u/v/weiv_salttr3d : the components of the salt transport due to the Gent-Mc Williams? scheme.

File size: 19.7 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                    ! store now fields before applying the Asselin filter
130         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
131         ! Asselin filter trend
132         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
133         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
134         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
135            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
136            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
137         ENDIF
138         ! total trend (note slightly inaccurate because tsb is time-filtered and tsa isn't)
139         ztrdt(:,:,:) = 0._wp
140         ztrds(:,:,:) = 0._wp
141         DO jk = 1, jpkm1
142            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / r2dtra(jk) 
143            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) 
144         END DO
145         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
146         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
147      ENDIF
148
149      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
150         DO jn = 1, jpts
151            DO jk = 1, jpkm1
152               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
153            END DO
154         END DO
155      ELSE                                            ! Leap-Frog + Asselin filter time stepping
156         !
157         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   &
158           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)
159         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
160         ENDIF
161      ENDIF     
162      !
163     ! trends computation
164      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
165         DO jk = 1, jpkm1
166            zfact = 1._wp / r2dtra(jk)             
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      !!              - save in (ta,sa) an average over the three time levels
193      !!             which will be used to compute rdn and thus the semi-implicit
194      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
195      !!              - swap tracer fields to prepare the next time_step.
196      !!                This can be summurized for tempearture as:
197      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
198      !!             ztm = 0                                   otherwise
199      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
200      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
201      !!             tn  = ta 
202      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
203      !!
204      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
205      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
206      !!----------------------------------------------------------------------
207      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
208      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
209      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
210      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
211      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
212      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
213      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
214      !
215      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
216      LOGICAL  ::   ll_tra_hpg       ! local logical
217      REAL(wp) ::   ztn, ztd         ! local scalars
218      !!----------------------------------------------------------------------
219
220      IF( kt == kit000 )  THEN
221         IF(lwp) WRITE(numout,*)
222         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
223         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
224      ENDIF
225      !
226      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
227      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
228      ENDIF
229      !
230      DO jn = 1, kjpt
231         !
232         DO jk = 1, jpkm1
233            DO jj = 1, jpj
234               DO ji = 1, jpi
235                  ztn = ptn(ji,jj,jk,jn)                                   
236                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
237                  !
238                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
239                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
240                  !
241                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
242               END DO
243           END DO
244         END DO
245         !
246      END DO
247      !
248   END SUBROUTINE tra_nxt_fix
249
250
251   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
252      !!----------------------------------------------------------------------
253      !!                   ***  ROUTINE tra_nxt_vvl  ***
254      !!
255      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
256      !!                and swap the tracer fields.
257      !!
258      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
259      !!              - save in (ta,sa) a thickness weighted average over the three
260      !!             time levels which will be used to compute rdn and thus the semi-
261      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
262      !!              - swap tracer fields to prepare the next time_step.
263      !!                This can be summurized for tempearture as:
264      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
265      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
266      !!             ztm = 0                                                       otherwise
267      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
268      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
269      !!             tn  = ta
270      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
271      !!
272      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
273      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
274      !!----------------------------------------------------------------------
275      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
276      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
277      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
278      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
279      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
280      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
281      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
282      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
283      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
284      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
285
286      !!     
287      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf, ll_isf   ! local logical
288      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
289      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
290      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
291      !!----------------------------------------------------------------------
292      !
293      IF( kt == kit000 )  THEN
294         IF(lwp) WRITE(numout,*)
295         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
296         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
297      ENDIF
298      !
299      IF( cdtype == 'TRA' )  THEN   
300         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
301         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
302         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
303         IF (nn_isf .GE. 1) THEN
304            ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing
305         ELSE
306            ll_isf = .FALSE.
307         END IF
308      ELSE                         
309         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
310         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
311         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
312         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing
313      ENDIF
314      !
315      DO jn = 1, kjpt     
316         DO jk = 1, jpkm1
317            zfact1 = atfp * p2dt(jk)
318            zfact2 = zfact1 / rau0
319            DO jj = 1, jpj
320               DO ji = 1, jpi
321                  ze3t_b = fse3t_b(ji,jj,jk)
322                  ze3t_n = fse3t_n(ji,jj,jk)
323                  ze3t_a = fse3t_a(ji,jj,jk)
324                  !                                         ! tracer content at Before, now and after
325                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
326                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
327                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
328                  !
329                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
330                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
331                  !
332                  ze3t_f = ze3t_n + atfp * ze3t_d
333                  ztc_f  = ztc_n  + atfp * ztc_d
334                  !
335                  IF( jk == mikt(ji,jj) ) THEN           ! first level
336                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
337                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
338                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
339                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
340                  ENDIF
341
342                  ! solar penetration (temperature only)
343                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
344                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
345
346                  ! river runoff
347                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
348                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
349                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
350
351                  ! ice shelf
352                  IF( ll_isf ) THEN
353                     ! level fully include in the Losch_2008 ice shelf boundary layer
354                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
355                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
356                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
357                     ! level partially include in Losch_2008 ice shelf boundary layer
358                     IF ( jk == misfkb(ji,jj) )                                                   &
359                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
360                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
361                  END IF
362
363                  ze3t_f = 1.e0 / ze3t_f
364                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
365                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
366                  !
367                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
368                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
369                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
370                  ENDIF
371               END DO
372            END DO
373         END DO
374         !
375      END DO
376      !
377   END SUBROUTINE tra_nxt_vvl
378
379   !!======================================================================
380END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.