source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 4328

Last change on this file since 4328 was 4328, checked in by davestorkey, 7 years ago

Remove OBC module at NEMO 3.6. See ticket #1189.

  • Property svn:keywords set to Id
File size: 16.9 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 zdf_oce         ! ???
30   USE domvvl          ! variable volume
31   USE dynspg_oce      ! surface     pressure gradient variables
32   USE dynhpg          ! hydrostatic pressure gradient
33   USE trdmod_oce      ! ocean space and time domain variables
34   USE trdtra          ! ocean active tracers trends
35   USE phycst
36   USE bdy_oce
37   USE bdytra          ! open boundary condition (bdy_tra routine)
38   USE in_out_manager  ! I/O manager
39   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
40   USE prtctl          ! Print control
41   USE traqsr          ! penetrative solar radiation (needed for nksr)
42#if defined key_agrif
43   USE agrif_opa_update
44   USE agrif_opa_interp
45#endif
46   USE wrk_nemo        ! Memory allocation
47   USE timing          ! Timing
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC   tra_nxt       ! routine called by step.F90
53   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
54   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
55
56   REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE tra_nxt( kt )
68      !!----------------------------------------------------------------------
69      !!                   ***  ROUTINE tranxt  ***
70      !!
71      !! ** Purpose :   Apply the boundary condition on the after temperature 
72      !!             and salinity fields, achieved the time stepping by adding
73      !!             the Asselin filter on now fields and swapping the fields.
74      !!
75      !! ** Method  :   At this stage of the computation, ta and sa are the
76      !!             after temperature and salinity as the time stepping has
77      !!             been performed in trazdf_imp or trazdf_exp module.
78      !!
79      !!              - Apply lateral boundary conditions on (ta,sa)
80      !!             at the local domain   boundaries through lbc_lnk call,
81      !!             at the one-way open boundaries (lk_bdy=T),
82      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
83      !!
84      !!              - Update lateral boundary conditions on AGRIF children
85      !!             domains (lk_agrif=T)
86      !!
87      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
88      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
89      !!----------------------------------------------------------------------
90      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
91      !!
92      INTEGER  ::   jk, jn    ! dummy loop indices
93      REAL(wp) ::   zfact     ! local scalars
94      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
95      !!----------------------------------------------------------------------
96      !
97      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt')
98      !
99      IF( kt == nit000 ) THEN
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
102         IF(lwp) WRITE(numout,*) '~~~~~~~'
103         !
104         rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp)      ! Brown & Campana parameter for semi-implicit hpg
105      ENDIF
106
107      ! Update after tracer on domain lateral boundaries
108      !
109      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign)
110      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
111      !
112#if defined key_bdy 
113      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
114#endif
115#if defined key_agrif
116      CALL Agrif_tra                     ! AGRIF zoom boundaries
117#endif
118 
119      ! set time step size (Euler/Leapfrog)
120      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
121      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
122      ENDIF
123
124      ! trends computation initialisation
125      IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter
126         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
127         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
128         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
129      ENDIF
130
131      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
132         DO jn = 1, jpts
133            DO jk = 1, jpkm1
134               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
135            END DO
136         END DO
137      ELSE                                            ! Leap-Frog + Asselin filter time stepping
138         !
139         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
140         ELSE                 ;   CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
141         ENDIF
142      ENDIF 
143      !
144#if defined key_agrif
145      ! Update tracer at AGRIF zoom boundaries
146      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
147#endif     
148      !
149      ! trends computation
150      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
151         DO jk = 1, jpkm1
152            zfact = 1.e0_wp / r2dtra(jk)             
153            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
154            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
155         END DO
156         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt )
157         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds )
158         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
159      END IF
160      !
161      !                        ! control print
162      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
163         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
164      !
165      !
166      IF( nn_timing == 1 )  CALL timing_stop('tra_nxt')
167      !
168   END SUBROUTINE tra_nxt
169
170
171   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
172      !!----------------------------------------------------------------------
173      !!                   ***  ROUTINE tra_nxt_fix  ***
174      !!
175      !! ** Purpose :   fixed volume: apply the Asselin time filter and
176      !!                swap the tracer fields.
177      !!
178      !! ** Method  : - Apply a Asselin time filter on now fields.
179      !!              - save in (ta,sa) an average over the three time levels
180      !!             which will be used to compute rdn and thus the semi-implicit
181      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
182      !!              - swap tracer fields to prepare the next time_step.
183      !!                This can be summurized for tempearture as:
184      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
185      !!             ztm = 0                                   otherwise
186      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
187      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
188      !!             tn  = ta 
189      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
190      !!
191      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
192      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
193      !!----------------------------------------------------------------------
194      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
195      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
196      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
197      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
198      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
199      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
200      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
201      !
202      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
203      LOGICAL  ::   ll_tra_hpg       ! local logical
204      REAL(wp) ::   ztn, ztd         ! local scalars
205      !!----------------------------------------------------------------------
206
207      IF( kt == kit000 )  THEN
208         IF(lwp) WRITE(numout,*)
209         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
210         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
211      ENDIF
212      !
213      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
214      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
215      ENDIF
216      !
217      DO jn = 1, kjpt
218         !
219         DO jk = 1, jpkm1
220            DO jj = 1, jpj
221               DO ji = 1, jpi
222                  ztn = ptn(ji,jj,jk,jn)                                   
223                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
224                  !
225                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
226                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
227                  !
228                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
229               END DO
230           END DO
231         END DO
232         !
233      END DO
234      !
235   END SUBROUTINE tra_nxt_fix
236
237
238   SUBROUTINE tra_nxt_vvl( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
239      !!----------------------------------------------------------------------
240      !!                   ***  ROUTINE tra_nxt_vvl  ***
241      !!
242      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
243      !!                and swap the tracer fields.
244      !!
245      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
246      !!              - save in (ta,sa) a thickness weighted average over the three
247      !!             time levels which will be used to compute rdn and thus the semi-
248      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
249      !!              - swap tracer fields to prepare the next time_step.
250      !!                This can be summurized for tempearture as:
251      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
252      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
253      !!             ztm = 0                                                       otherwise
254      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
255      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
256      !!             tn  = ta
257      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
258      !!
259      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
260      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
261      !!----------------------------------------------------------------------
262      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
263      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
264      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
265      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
266      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
267      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
268      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
269      !!     
270      LOGICAL  ::   ll_tra, ll_tra_hpg, ll_traqsr   ! local logical
271      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
272      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
273      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
274      !!----------------------------------------------------------------------
275      !
276      IF( kt == kit000 )  THEN
277         IF(lwp) WRITE(numout,*)
278         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
279         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
280      ENDIF
281      !
282      IF( cdtype == 'TRA' )  THEN   
283         ll_tra     = .TRUE.           ! active tracers case 
284         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
285         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
286      ELSE                         
287         ll_tra     = .FALSE.          ! passive tracers case
288         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
289         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
290      ENDIF
291      !
292      DO jn = 1, kjpt     
293         DO jk = 1, jpkm1
294            zfact1 = atfp * rdttra(jk)
295            zfact2 = zfact1 / rau0
296            DO jj = 1, jpj
297               DO ji = 1, jpi
298                  ze3t_b = fse3t_b(ji,jj,jk)
299                  ze3t_n = fse3t_n(ji,jj,jk)
300                  ze3t_a = fse3t_a(ji,jj,jk)
301                  !                                         ! tracer content at Before, now and after
302                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
303                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
304                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
305                  !
306                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
307                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
308                  !
309                  ze3t_f = ze3t_n + atfp * ze3t_d
310                  ztc_f  = ztc_n  + atfp * ztc_d
311                  !
312                  IF( ll_tra .AND. jk == 1 ) THEN           ! first level only for T & S
313                      ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
314                      ztc_f  = ztc_f  - zfact1 * ( sbc_tsc(ji,jj,jn) - sbc_tsc_b(ji,jj,jn) )
315                  ENDIF
316                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
317                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
318
319                   ze3t_f = 1.e0 / ze3t_f
320                   ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
321                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
322                   !
323                   IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
324                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
325                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
326                   ENDIF
327               END DO
328            END DO
329         END DO
330         !
331      END DO
332      !
333   END SUBROUTINE tra_nxt_vvl
334
335   !!======================================================================
336END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.