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

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 7256

Last change on this file since 7256 was 7256, checked in by cbricaud, 7 years ago

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