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

source: branches/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 10207

Last change on this file since 10207 was 10207, checked in by cmao, 5 years ago

remove svn keyword
File size: 28.3 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   USE crs
52
53   IMPLICIT NONE
54   PRIVATE
55
56   PUBLIC   tra_nxt       ! routine called by step.F90
57   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
58   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
59   PUBLIC   tra_nxt_vvl_crs ! to be used in trcnxt
60
61   REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg
62
63   !! * Substitutions
64#  include "domzgr_substitute.h90"
65   !!----------------------------------------------------------------------
66   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
67   !! $Id$
68   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
69   !!----------------------------------------------------------------------
70CONTAINS
71
72   SUBROUTINE tra_nxt( kt )
73      !!----------------------------------------------------------------------
74      !!                   ***  ROUTINE tranxt  ***
75      !!
76      !! ** Purpose :   Apply the boundary condition on the after temperature 
77      !!             and salinity fields, achieved the time stepping by adding
78      !!             the Asselin filter on now fields and swapping the fields.
79      !!
80      !! ** Method  :   At this stage of the computation, ta and sa are the
81      !!             after temperature and salinity as the time stepping has
82      !!             been performed in trazdf_imp or trazdf_exp module.
83      !!
84      !!              - Apply lateral boundary conditions on (ta,sa)
85      !!             at the local domain   boundaries through lbc_lnk call,
86      !!             at the one-way open boundaries (lk_bdy=T),
87      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
88      !!
89      !!              - Update lateral boundary conditions on AGRIF children
90      !!             domains (lk_agrif=T)
91      !!
92      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
93      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
94      !!----------------------------------------------------------------------
95      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
96      !!
97      INTEGER  ::   jk, jn    ! dummy loop indices
98      REAL(wp) ::   zfact     ! local scalars
99      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
100      !!----------------------------------------------------------------------
101      !
102      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt')
103      !
104      IF( kt == nit000 ) THEN
105         IF(lwp) WRITE(numout,*)
106         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
107         IF(lwp) WRITE(numout,*) '~~~~~~~'
108         !
109         rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp)      ! Brown & Campana parameter for semi-implicit hpg
110      ENDIF
111
112      ! Update after tracer on domain lateral boundaries
113      !
114#if defined key_agrif
115      CALL Agrif_tra                     ! AGRIF zoom boundaries
116#endif
117      !
118      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign)
119      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
120      !
121#if defined key_bdy 
122      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
123#endif
124 
125      ! set time step size (Euler/Leapfrog)
126      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
127      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
128      ENDIF
129
130      ! trends computation initialisation
131      IF( l_trdtra )   THEN                   
132         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
133         ztrdt(:,:,jpk) = 0._wp
134         ztrds(:,:,jpk) = 0._wp
135         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
136            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
137            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
138         ENDIF
139         ! total trend for the non-time-filtered variables.
140         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms
141         IF( lk_vvl ) THEN
142            DO jk = 1, jpkm1
143               zfact = 1.0 / rdttra(jk)
144               ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - tsn(:,:,jk,jp_tem)) * zfact
145               ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - tsn(:,:,jk,jp_sal)) * zfact
146            END DO
147         ELSE
148            DO jk = 1, jpkm1
149               zfact = 1.0 / rdttra(jk)
150               ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsn(:,:,jk,jp_tem) ) * zfact 
151               ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsn(:,:,jk,jp_sal) ) * zfact 
152            END DO
153         END IF
154         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
155         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
156         IF( .NOT.lk_vvl )  THEN
157            ! Store now fields before applying the Asselin filter
158            ! in order to calculate Asselin filter trend later.
159            ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
160            ztrds(:,:,:) = tsn(:,:,:,jp_sal)
161         END IF
162      ENDIF
163
164      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
165         DO jn = 1, jpts
166            DO jk = 1, jpkm1
167               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
168            END DO
169         END DO
170         IF (l_trdtra.AND.lk_vvl) THEN      ! Zero Asselin filter contribution must be explicitly written out since for vvl
171                                            ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
172            ztrdt(:,:,:) = 0._wp
173            ztrds(:,:,:) = 0._wp
174            CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
175            CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
176         END IF
177      ELSE                                            ! Leap-Frog + Asselin filter time stepping
178         !
179         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   &
180           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)
181         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
182         ENDIF
183      ENDIF     
184      !
185     ! trends computation
186      IF( l_trdtra.AND..NOT.lk_vvl) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
187         DO jk = 1, jpkm1
188            zfact = 1._wp / r2dtra(jk)             
189            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
190            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
191         END DO
192         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
193         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
194      END IF
195      IF( l_trdtra) CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
196      !
197      !                        ! control print
198      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
199         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
200      !
201      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
202      !
203   END SUBROUTINE tra_nxt
204
205
206   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
207      !!----------------------------------------------------------------------
208      !!                   ***  ROUTINE tra_nxt_fix  ***
209      !!
210      !! ** Purpose :   fixed volume: apply the Asselin time filter and
211      !!                swap the tracer fields.
212      !!
213      !! ** Method  : - Apply a Asselin time filter on now fields.
214      !!              - save in (ta,sa) an average over the three time levels
215      !!             which will be used to compute rdn and thus the semi-implicit
216      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
217      !!              - swap tracer fields to prepare the next time_step.
218      !!                This can be summurized for tempearture as:
219      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
220      !!             ztm = 0                                   otherwise
221      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
222      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
223      !!             tn  = ta 
224      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
225      !!
226      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
227      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
228      !!----------------------------------------------------------------------
229      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
230      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
231      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
232      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
233      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
234      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
235      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
236      !
237      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
238      LOGICAL  ::   ll_tra_hpg       ! local logical
239      REAL(wp) ::   ztn, ztd         ! local scalars
240      !!----------------------------------------------------------------------
241
242      IF( kt == kit000 )  THEN
243         IF(lwp) WRITE(numout,*)
244         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
245         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
246      ENDIF
247      !
248      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
249      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
250      ENDIF
251      !
252      DO jn = 1, kjpt
253         !
254         DO jk = 1, jpkm1
255            DO jj = 1, jpj
256               DO ji = 1, jpi
257                  ztn = ptn(ji,jj,jk,jn)                                   
258                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
259                  !
260                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
261                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
262                  !
263                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
264               END DO
265           END DO
266         END DO
267         !
268      END DO
269      !
270   END SUBROUTINE tra_nxt_fix
271
272
273   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
274      !!----------------------------------------------------------------------
275      !!                   ***  ROUTINE tra_nxt_vvl  ***
276      !!
277      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
278      !!                and swap the tracer fields.
279      !!
280      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
281      !!              - save in (ta,sa) a thickness weighted average over the three
282      !!             time levels which will be used to compute rdn and thus the semi-
283      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
284      !!              - swap tracer fields to prepare the next time_step.
285      !!                This can be summurized for tempearture as:
286      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
287      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
288      !!             ztm = 0                                                       otherwise
289      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
290      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
291      !!             tn  = ta
292      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
293      !!
294      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
295      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
296      !!----------------------------------------------------------------------
297      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
298      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
299      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
300      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
301      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
302      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
303      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
304      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
305      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
306      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
307
308      !!     
309      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf, ll_isf   ! local logical
310      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
311      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
312      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
313      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrd_atf
314      !!----------------------------------------------------------------------
315      !
316      IF( kt == kit000 )  THEN
317         IF(lwp) WRITE(numout,*)
318         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
319         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
320      ENDIF
321      !
322      IF( cdtype == 'TRA' )  THEN   
323         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
324         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
325         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
326         IF (nn_isf .GE. 1) THEN
327            ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing
328         ELSE
329            ll_isf = .FALSE.
330         END IF
331      ELSE                         
332         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
333         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
334         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
335         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing
336      ENDIF
337      !
338      IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) )   THEN
339         CALL wrk_alloc( jpi, jpj, jpk, kjpt, ztrd_atf )
340         ztrd_atf(:,:,:,:) = 0.0_wp
341      ENDIF
342      DO jn = 1, kjpt     
343         DO jk = 1, jpkm1
344            zfact = 1._wp / r2dtra(jk)
345            zfact1 = atfp * p2dt(jk)
346            zfact2 = zfact1 / rau0
347            DO jj = 1, jpj
348               DO ji = 1, jpi
349                  ze3t_b = fse3t_b(ji,jj,jk)
350                  ze3t_n = fse3t_n(ji,jj,jk)
351                  ze3t_a = fse3t_a(ji,jj,jk)
352                  !                                         ! tracer content at Before, now and after
353                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
354                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
355                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
356                  !
357                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
358                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
359                  !
360                  ze3t_f = ze3t_n + atfp * ze3t_d
361                  ztc_f  = ztc_n  + atfp * ztc_d
362                  !
363                  IF( jk == mikt(ji,jj) ) THEN           ! first level
364                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
365                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
366                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
367                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
368                  ENDIF
369
370                  ! solar penetration (temperature only)
371                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
372                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
373
374                  ! river runoff
375                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
376                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
377                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
378
379                  ! ice shelf
380                  IF( ll_isf ) THEN
381                     ! level fully include in the Losch_2008 ice shelf boundary layer
382                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
383                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
384                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
385                     ! level partially include in Losch_2008 ice shelf boundary layer
386                     IF ( jk == misfkb(ji,jj) )                                                   &
387                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
388                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
389                  END IF
390
391                  ze3t_f = 1.e0 / ze3t_f
392                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
393                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
394                  !
395                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
396                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
397                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
398                  ENDIF
399                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
400                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
401                  ENDIF
402               END DO
403            END DO
404         END DO
405         !
406      END DO
407      !
408      IF( l_trdtra .and. cdtype == 'TRA' ) THEN
409         CALL trd_tra( kt, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
410         CALL trd_tra( kt, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
411         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
412      ENDIF
413      IF( l_trdtrc .and. cdtype == 'TRC' ) THEN
414         DO jn = 1, kjpt
415            CALL trd_tra( kt, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
416         END DO
417         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
418      ENDIF
419
420   END SUBROUTINE tra_nxt_vvl
421
422  SUBROUTINE tra_nxt_vvl_crs( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
423      !!----------------------------------------------------------------------
424      !!                   ***  ROUTINE tra_nxt_vvl  ***
425      !!
426      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
427      !!                and swap the tracer fields.
428      !!
429      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
430      !!              - save in (ta,sa) a thickness weighted average over the three
431      !!             time levels which will be used to compute rdn and thus the semi-
432      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
433      !!              - swap tracer fields to prepare the next time_step.
434      !!                This can be summurized for tempearture as:
435      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
436      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
437      !!             ztm = 0                                                       otherwise
438      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
439      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
440      !!             tn  = ta
441      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
442      !!
443      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
444      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
445      !!----------------------------------------------------------------------
446      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
447      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
448      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
449      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
450      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
451      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
452      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
453      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
454      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
455      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
456
457      !!     
458      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical
459      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
460      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
461      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
462      !!----------------------------------------------------------------------
463      !!----------------------------------------------------------------------
464      !
465      IF( kt == kit000 )  THEN
466         IF(lwp) WRITE(numout,*)
467         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
468         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
469      ENDIF
470      !
471      IF( cdtype == 'TRA' )  THEN
472         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
473         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
474         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
475      ELSE
476         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
477         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
478         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
479      ENDIF
480      !
481      DO jn = 1, kjpt
482         DO jk = 1, jpkm1
483            zfact1 = atfp * p2dt(jk)
484            zfact2 = zfact1 / rau0
485            DO jj = 1, jpj
486               DO ji = 1, jpi
487                  ze3t_b = fse3t_b_crs(ji,jj,jk)
488                  ze3t_n = fse3t_n_crs(ji,jj,jk)
489                  ze3t_a = fse3t_a_crs(ji,jj,jk)
490                  !                                         ! tracer content at Before, now and after
491                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
492                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
493                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
494                  !
495                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
496                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
497                  !
498                  ze3t_f = ze3t_n + atfp * ze3t_d
499                  ztc_f  = ztc_n  + atfp * ztc_d
500                  !
501                  IF( jk == mikt_crs(ji,jj) ) THEN           ! first level
502                     ze3t_f = ze3t_f - zfact2 * ( (emp_b_crs(ji,jj)    - emp_crs(ji,jj)   )  &
503                            &                   - (rnf_b_crs(ji,jj)    - rnf_crs(ji,jj)   )  &
504                            &                   + (fwfisf_b_crs(ji,jj) - fwfisf_crs(ji,jj))  )
505                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
506                  ENDIF
507
508                  !cbr: not adapted for coarsening for the moment
509                  ! ice shelf
510                  !IF( ll_isf ) THEN
511                  !   ! level fully include in the Losch_2008 ice shelf boundary layer
512                  !   IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
513                  !      ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
514                  !             &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
515                  !   ! level partially include in Losch_2008 ice shelf boundary layer
516                  !   IF ( jk == misfkb(ji,jj) )                                                   &
517                  !      ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
518                  !             &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
519                  !END IF
520
521                  ze3t_f = 1.e0 / ze3t_f
522                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
523                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
524                  !
525               END DO
526            END DO
527         END DO
528         !
529      END DO
530      !
531   END SUBROUTINE tra_nxt_vvl_crs
532
533
534   !!======================================================================
535END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.