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 @ 6772

Last change on this file since 6772 was 6772, checked in by cbricaud, 8 years ago

clean in coarsening branch
  • Property svn:keywords set to Id
File size: 24.2 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 zdf_oce         ! ocean vertical mixing
31   USE domvvl          ! variable volume
32   USE dynspg_oce      ! surface     pressure gradient variables
33   USE dynhpg          ! hydrostatic pressure gradient
34   USE trd_oce         ! trends: ocean variables
35   USE trdtra          ! trends manager: tracers
36   USE traqsr          ! penetrative solar radiation (needed for nksr)
37   USE phycst          ! physical constant
38   USE ldftra_oce      ! lateral physics on tracers
39   USE bdy_oce         ! BDY open boundary condition variables
40   USE bdytra          ! open boundary condition (bdy_tra routine)
41   !
42   USE in_out_manager  ! I/O manager
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE prtctl          ! Print control
45   USE wrk_nemo        ! Memory allocation
46   USE timing          ! Timing
47#if defined key_agrif
48   USE agrif_opa_update
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      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign)
115      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
116      !
117#if defined key_bdy 
118      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
119#endif
120#if defined key_agrif
121      CALL Agrif_tra                     ! AGRIF zoom boundaries
122#endif
123 
124      ! set time step size (Euler/Leapfrog)
125      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
126      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
127      ENDIF
128
129      ! trends computation initialisation
130      IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter
131         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
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      ENDIF
139
140      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
141         DO jn = 1, jpts
142            DO jk = 1, jpkm1
143               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
144            END DO
145         END DO
146      ELSE                                            ! Leap-Frog + Asselin filter time stepping
147         !
148         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   &
149           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)
150         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
151         ENDIF
152      ENDIF 
153      !
154#if defined key_agrif
155      ! Update tracer at AGRIF zoom boundaries
156      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
157#endif     
158      !
159      ! trends computation
160      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
161         DO jk = 1, jpkm1
162            zfact = 1._wp / r2dtra(jk)             
163            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
164            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
165         END DO
166         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
167         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
168         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
169      END IF
170      !
171      !                        ! control print
172      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
173         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
174      !
175      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
176      !
177   END SUBROUTINE tra_nxt
178
179
180   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
181      !!----------------------------------------------------------------------
182      !!                   ***  ROUTINE tra_nxt_fix  ***
183      !!
184      !! ** Purpose :   fixed volume: apply the Asselin time filter and
185      !!                swap the tracer fields.
186      !!
187      !! ** Method  : - Apply a Asselin time filter on now fields.
188      !!              - save in (ta,sa) an average over the three time levels
189      !!             which will be used to compute rdn and thus the semi-implicit
190      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
191      !!              - swap tracer fields to prepare the next time_step.
192      !!                This can be summurized for tempearture as:
193      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
194      !!             ztm = 0                                   otherwise
195      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
196      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
197      !!             tn  = ta 
198      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
199      !!
200      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
201      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
202      !!----------------------------------------------------------------------
203      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
204      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
205      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
206      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
207      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
208      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
209      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
210      !
211      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
212      LOGICAL  ::   ll_tra_hpg       ! local logical
213      REAL(wp) ::   ztn, ztd         ! local scalars
214      !!----------------------------------------------------------------------
215
216      IF( kt == kit000 )  THEN
217         IF(lwp) WRITE(numout,*)
218         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
219         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
220      ENDIF
221      !
222      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
223      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
224      ENDIF
225      !
226      DO jn = 1, kjpt
227         !
228         DO jk = 1, jpkm1
229            DO jj = 1, jpj
230               DO ji = 1, jpi
231                  ztn = ptn(ji,jj,jk,jn)                                   
232                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
233                  !
234                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
235                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
236                  !
237                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
238               END DO
239           END DO
240         END DO
241         !
242      END DO
243      !
244   END SUBROUTINE tra_nxt_fix
245
246
247   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
248      !!----------------------------------------------------------------------
249      !!                   ***  ROUTINE tra_nxt_vvl  ***
250      !!
251      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
252      !!                and swap the tracer fields.
253      !!
254      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
255      !!              - save in (ta,sa) a thickness weighted average over the three
256      !!             time levels which will be used to compute rdn and thus the semi-
257      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
258      !!              - swap tracer fields to prepare the next time_step.
259      !!                This can be summurized for tempearture as:
260      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
261      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
262      !!             ztm = 0                                                       otherwise
263      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
264      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
265      !!             tn  = ta
266      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
267      !!
268      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
269      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
270      !!----------------------------------------------------------------------
271      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
272      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
273      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
274      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
275      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
276      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
277      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
278      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
279      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
280      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
281
282      !!     
283      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical
284      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
285      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
286      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
287      !!----------------------------------------------------------------------
288      !
289      IF( kt == kit000 )  THEN
290         IF(lwp) WRITE(numout,*)
291         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
292         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
293      ENDIF
294      !
295      IF( cdtype == 'TRA' )  THEN   
296         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
297         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
298         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
299      ELSE                         
300         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
301         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
302         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
303      ENDIF
304      !
305      DO jn = 1, kjpt     
306         DO jk = 1, jpkm1
307            zfact1 = atfp * p2dt(jk)
308            zfact2 = zfact1 / rau0
309            DO jj = 1, jpj
310               DO ji = 1, jpi
311                  ze3t_b = fse3t_b(ji,jj,jk)
312                  ze3t_n = fse3t_n(ji,jj,jk)
313                  ze3t_a = fse3t_a(ji,jj,jk)
314                  !                                         ! tracer content at Before, now and after
315                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
316                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
317                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
318                  !
319                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
320                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
321                  !
322                  ze3t_f = ze3t_n + atfp * ze3t_d
323                  ztc_f  = ztc_n  + atfp * ztc_d
324                  !
325                  IF( jk == 1 ) THEN           ! first level
326                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) )
327                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
328                  ENDIF
329
330                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
331                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
332
333                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )   &            ! river runoffs
334                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
335                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
336
337                  ze3t_f = 1.e0 / ze3t_f
338                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
339                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
340                  !
341                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
342                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
343                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
344                  ENDIF
345               END DO
346            END DO
347         END DO
348         !
349      END DO
350      !
351   END SUBROUTINE tra_nxt_vvl
352
353  SUBROUTINE tra_nxt_vvl_crs( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
354      !!----------------------------------------------------------------------
355      !!                   ***  ROUTINE tra_nxt_vvl  ***
356      !!
357      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
358      !!                and swap the tracer fields.
359      !!
360      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
361      !!              - save in (ta,sa) a thickness weighted average over the three
362      !!             time levels which will be used to compute rdn and thus the semi-
363      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
364      !!              - swap tracer fields to prepare the next time_step.
365      !!                This can be summurized for tempearture as:
366      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
367      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
368      !!             ztm = 0                                                       otherwise
369      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
370      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
371      !!             tn  = ta
372      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
373      !!
374      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
375      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
376      !!----------------------------------------------------------------------
377      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
378      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
379      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
380      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
381      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
382      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
383      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
384      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
385      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
386      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
387
388      !!     
389      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical
390      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
391      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
392      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
393      !!----------------------------------------------------------------------
394      !!----------------------------------------------------------------------
395      !
396      IF( kt == kit000 )  THEN
397         IF(lwp) WRITE(numout,*)
398         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
399         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
400      ENDIF
401      !
402      IF( cdtype == 'TRA' )  THEN
403         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
404         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
405         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
406      ELSE
407         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
408         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
409         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
410      ENDIF
411      !
412      DO jn = 1, kjpt
413         DO jk = 1, jpkm1
414            zfact1 = atfp * p2dt(jk)
415            zfact2 = zfact1 / rau0
416            DO jj = 1, jpj
417               DO ji = 1, jpi
418                  ze3t_b = fse3t_b_crs(ji,jj,jk)
419                  ze3t_n = fse3t_n_crs(ji,jj,jk)
420                  ze3t_a = fse3t_a_crs(ji,jj,jk)
421                  !                                         ! tracer content at Before, now and after
422                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
423                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
424                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
425                  !
426                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
427                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
428                  !
429                  ze3t_f = ze3t_n + atfp * ze3t_d
430                  ztc_f  = ztc_n  + atfp * ztc_d
431                  !
432                  IF( jk == 1 ) THEN           ! first level
433                     ze3t_f = ze3t_f - zfact2 * ( emp_b_crs(ji,jj) - emp_crs(ji,jj) + rnf_crs(ji,jj) - rnf_b_crs(ji,jj) )
434                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
435                  ENDIF
436!cbr as it is a subroutine dedicated to crs, TRA options are not necessary
437!cbr                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
438!cbr                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )
439!cbr
440!cbr                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )   &                  ! river runoffs
441!cbr                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) &
442!cbr                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
443
444                  ze3t_f = 1.e0 / ze3t_f
445                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
446                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
447                  !
448                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
449                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
450                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
451                  ENDIF
452               END DO
453            END DO
454         END DO
455         !
456      END DO
457      !
458   END SUBROUTINE tra_nxt_vvl_crs
459
460
461   !!======================================================================
462END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.