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

source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 8485

Last change on this file since 8485 was 8333, checked in by timgraham, 7 years ago

GMED ticket #336 - fix for bug in Asselin filter trend when TOP is enabled

File size: 21.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(:,:,jpk) = 0._wp
132         ztrds(:,:,jpk) = 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         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms
139         IF( lk_vvl ) THEN
140            DO jk = 1, jpkm1
141               zfact = 1.0 / rdttra(jk)
142               ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - tsn(:,:,jk,jp_tem)) * zfact
143               ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - tsn(:,:,jk,jp_sal)) * zfact
144            END DO
145         ELSE
146            DO jk = 1, jpkm1
147               zfact = 1.0 / rdttra(jk)
148               ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsn(:,:,jk,jp_tem) ) * zfact 
149               ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsn(:,:,jk,jp_sal) ) * zfact 
150            END DO
151         END IF
152         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
153         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
154         IF( .NOT.lk_vvl )  THEN
155            ! Store now fields before applying the Asselin filter
156            ! in order to calculate Asselin filter trend later.
157            ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
158            ztrds(:,:,:) = tsn(:,:,:,jp_sal)
159         END IF
160      ENDIF
161
162      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
163         DO jn = 1, jpts
164            DO jk = 1, jpkm1
165               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
166            END DO
167         END DO
168         IF (l_trdtra.AND.lk_vvl) THEN      ! Zero Asselin filter contribution must be explicitly written out since for vvl
169                                            ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
170            ztrdt(:,:,:) = 0._wp
171            ztrds(:,:,:) = 0._wp
172            CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
173            CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
174         END IF
175      ELSE                                            ! Leap-Frog + Asselin filter time stepping
176         !
177         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   &
178           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)
179         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
180         ENDIF
181      ENDIF     
182      !
183     ! trends computation
184      IF( l_trdtra.AND..NOT.lk_vvl) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
185         DO jk = 1, jpkm1
186            zfact = 1._wp / r2dtra(jk)             
187            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
188            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
189         END DO
190         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
191         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
192      END IF
193      IF( l_trdtra) CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
194      !
195      !                        ! control print
196      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
197         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
198      !
199      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
200      !
201   END SUBROUTINE tra_nxt
202
203
204   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
205      !!----------------------------------------------------------------------
206      !!                   ***  ROUTINE tra_nxt_fix  ***
207      !!
208      !! ** Purpose :   fixed volume: apply the Asselin time filter and
209      !!                swap the tracer fields.
210      !!
211      !! ** Method  : - Apply a Asselin time filter on now fields.
212      !!              - save in (ta,sa) an average over the three time levels
213      !!             which will be used to compute rdn and thus the semi-implicit
214      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
215      !!              - swap tracer fields to prepare the next time_step.
216      !!                This can be summurized for tempearture as:
217      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
218      !!             ztm = 0                                   otherwise
219      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
220      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
221      !!             tn  = ta 
222      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
223      !!
224      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
225      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
226      !!----------------------------------------------------------------------
227      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
228      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
229      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
230      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
231      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
232      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
233      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
234      !
235      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
236      LOGICAL  ::   ll_tra_hpg       ! local logical
237      REAL(wp) ::   ztn, ztd         ! local scalars
238      !!----------------------------------------------------------------------
239
240      IF( kt == kit000 )  THEN
241         IF(lwp) WRITE(numout,*)
242         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
243         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
244      ENDIF
245      !
246      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
247      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
248      ENDIF
249      !
250      DO jn = 1, kjpt
251         !
252         DO jk = 1, jpkm1
253            DO jj = 1, jpj
254               DO ji = 1, jpi
255                  ztn = ptn(ji,jj,jk,jn)                                   
256                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
257                  !
258                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
259                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
260                  !
261                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
262               END DO
263           END DO
264         END DO
265         !
266      END DO
267      !
268   END SUBROUTINE tra_nxt_fix
269
270
271   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
272      !!----------------------------------------------------------------------
273      !!                   ***  ROUTINE tra_nxt_vvl  ***
274      !!
275      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
276      !!                and swap the tracer fields.
277      !!
278      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
279      !!              - save in (ta,sa) a thickness weighted average over the three
280      !!             time levels which will be used to compute rdn and thus the semi-
281      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
282      !!              - swap tracer fields to prepare the next time_step.
283      !!                This can be summurized for tempearture as:
284      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
285      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
286      !!             ztm = 0                                                       otherwise
287      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
288      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
289      !!             tn  = ta
290      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
291      !!
292      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
293      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
294      !!----------------------------------------------------------------------
295      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
296      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
297      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
298      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
299      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
300      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
301      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
302      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
303      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
304      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
305
306      !!     
307      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf, ll_isf   ! local logical
308      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
309      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
310      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
311      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrd_atf
312      !!----------------------------------------------------------------------
313      !
314      IF( kt == kit000 )  THEN
315         IF(lwp) WRITE(numout,*)
316         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
317         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
318      ENDIF
319      !
320      IF( cdtype == 'TRA' )  THEN   
321         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
322         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
323         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
324         IF (nn_isf .GE. 1) THEN
325            ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing
326         ELSE
327            ll_isf = .FALSE.
328         END IF
329      ELSE                         
330         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
331         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
332         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
333         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing
334      ENDIF
335      !
336      IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) )   THEN
337         CALL wrk_alloc( jpi, jpj, jpk, kjpt, ztrd_atf )
338         ztrd_atf(:,:,:,:) = 0.0_wp
339      ENDIF
340      DO jn = 1, kjpt     
341         DO jk = 1, jpkm1
342            zfact = 1._wp / r2dtra(jk)
343            zfact1 = atfp * p2dt(jk)
344            zfact2 = zfact1 / rau0
345            DO jj = 1, jpj
346               DO ji = 1, jpi
347                  ze3t_b = fse3t_b(ji,jj,jk)
348                  ze3t_n = fse3t_n(ji,jj,jk)
349                  ze3t_a = fse3t_a(ji,jj,jk)
350                  !                                         ! tracer content at Before, now and after
351                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
352                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
353                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
354                  !
355                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
356                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
357                  !
358                  ze3t_f = ze3t_n + atfp * ze3t_d
359                  ztc_f  = ztc_n  + atfp * ztc_d
360                  !
361                  IF( jk == mikt(ji,jj) ) THEN           ! first level
362                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
363                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
364                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
365                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
366                  ENDIF
367
368                  ! solar penetration (temperature only)
369                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
370                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
371
372                  ! river runoff
373                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
374                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
375                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
376
377                  ! ice shelf
378                  IF( ll_isf ) THEN
379                     ! level fully include in the Losch_2008 ice shelf boundary layer
380                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
381                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
382                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
383                     ! level partially include in Losch_2008 ice shelf boundary layer
384                     IF ( jk == misfkb(ji,jj) )                                                   &
385                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
386                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
387                  END IF
388
389                  ze3t_f = 1.e0 / ze3t_f
390                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
391                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
392                  !
393                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
394                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
395                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
396                  ENDIF
397                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
398                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
399                  ENDIF
400               END DO
401            END DO
402         END DO
403         !
404      END DO
405      !
406      IF( l_trdtra .and. cdtype == 'TRA' ) THEN
407         CALL trd_tra( kt, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
408         CALL trd_tra( kt, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
409         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
410      ENDIF
411      IF( l_trdtrc .and. cdtype == 'TRC' ) THEN
412         DO jn = 1, kjpt
413            CALL trd_tra( kt, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
414         END DO
415         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
416      ENDIF
417
418   END SUBROUTINE tra_nxt_vvl
419
420   !!======================================================================
421END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.