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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranxt.F90 @ 2034

Last change on this file since 2034 was 2034, checked in by cetlod, 14 years ago

cosmetic changes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 20.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   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   tra_nxt       : time stepping on temperature and salinity
21   !!   tra_nxt_fix   : time stepping on temperature and salinity : fixed    volume case
22   !!   tra_nxt_vvl   : time stepping on temperature and salinity : variable volume case
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers variables
25   USE dom_oce         ! ocean space and time domain variables
26   USE zdf_oce         ! ???
27   USE domvvl          ! variable volume
28   USE dynspg_oce      ! surface     pressure gradient variables
29   USE dynhpg          ! hydrostatic pressure gradient
30   USE trdmod_oce      ! ocean space and time domain variables
31   USE trdtra          ! ocean active tracers trends
32   USE phycst
33   USE obctra          ! open boundary condition (obc_tra routine)
34   USE bdytra          ! Unstructured open boundary condition (bdy_tra routine)
35   USE in_out_manager  ! I/O manager
36   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
37   USE prtctl          ! Print control
38   USE traswp          ! swap array
39   USE agrif_opa_update
40   USE agrif_opa_interp
41   USE obc_oce
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   tra_nxt       ! routine called by step.F90
47   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
48   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
49
50   REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
51
52   !! * Substitutions
53#  include "domzgr_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
56   !! $Id$
57   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59
60CONTAINS
61
62   SUBROUTINE tra_nxt( kt )
63      !!----------------------------------------------------------------------
64      !!                   ***  ROUTINE tranxt  ***
65      !!
66      !! ** Purpose :   Apply the boundary condition on the after temperature 
67      !!             and salinity fields, achieved the time stepping by adding
68      !!             the Asselin filter on now fields and swapping the fields.
69      !!
70      !! ** Method  :   At this stage of the computation, ta and sa are the
71      !!             after temperature and salinity as the time stepping has
72      !!             been performed in trazdf_imp or trazdf_exp module.
73      !!
74      !!              - Apply lateral boundary conditions on (ta,sa)
75      !!             at the local domain   boundaries through lbc_lnk call,
76      !!             at the radiative open boundaries (lk_obc=T),
77      !!             at the relaxed   open boundaries (lk_bdy=T), and
78      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
79      !!
80      !!              - Update lateral boundary conditions on AGRIF children
81      !!             domains (lk_agrif=T)
82      !!
83      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
84      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
85      !!----------------------------------------------------------------------
86      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
87      !!
88      INTEGER  ::   jk    ! dummy loop indices
89      REAL(wp) ::   zfact ! temporary scalars
90      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
91
92      !!----------------------------------------------------------------------
93
94      IF( kt == nit000 ) THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
97         IF(lwp) WRITE(numout,*) '~~~~~~~'
98      ENDIF
99
100      ! Update after tracer on domain lateral boundaries
101      !
102      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )      ! local domain boundaries  (T-point, unchanged sign)
103      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
104      !
105#if defined key_obc || defined key_bdy || defined key_agrif
106      CALL tra_unswap
107#endif
108#if defined key_obc
109      IF( lk_obc )   CALL obc_tra( kt )  ! OBC open boundaries
110#endif
111#if defined key_bdy
112      CALL bdy_tra( kt )               ! BDY open boundaries
113#endif
114#if defined key_agrif
115      CALL Agrif_tra                   ! AGRIF zoom boundaries
116#endif
117#if defined key_obc || defined key_bdy || defined key_agrif
118      CALL tra_swap
119#endif
120 
121      ! set time step size (Euler/Leapfrog)
122      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
123      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
124      ENDIF
125
126      ! trends computation initialisation
127      IF( l_trdtra )   THEN                    !* store now fields before applying the Asselin filter
128         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
129         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsn(:,:,:,jp_sal)
130      ENDIF
131
132      ! Leap-Frog + Asselin filter time stepping
133      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, nit000, tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
134      ELSE                  ;   CALL tra_nxt_fix( kt, nit000, tsb, tsn, tsa, jpts )  ! fixed    volume level
135      ENDIF
136
137#if defined key_agrif
138      CALL tra_unswap
139      ! Update tracer at AGRIF zoom boundaries
140      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
141      CALL tra_swap
142#endif     
143
144      ! trends computation
145      IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt     
146         DO jk = 1, jpkm1
147            zfact = 1.e0 / r2dt_t(jk)             
148            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
149            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
150         END DO
151         CALL trd_tra( kt, 'TRA', jp_tem, jpt_trd_atf, ztrdt )
152         CALL trd_tra( kt, 'TRA', jp_sal, jpt_trd_atf, ztrds )
153         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
154      END IF
155
156      !                        ! control print
157      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
158         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
159      !
160   END SUBROUTINE tra_nxt
161
162   SUBROUTINE tra_nxt_fix( kt, kit000,                    &
163      &                               ptb, ptn, pta, kjpt  )
164      !!----------------------------------------------------------------------
165      !!                   ***  ROUTINE tra_nxt_fix  ***
166      !!
167      !! ** Purpose :   fixed volume: apply the Asselin time filter and
168      !!                swap the tracer fields.
169      !!
170      !! ** Method  : - Apply a Asselin time filter on now fields.
171      !!              - save in (ta,sa) an average over the three time levels
172      !!             which will be used to compute rdn and thus the semi-implicit
173      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
174      !!              - swap tracer fields to prepare the next time_step.
175      !!                This can be summurized for tempearture as:
176      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T
177      !!             ztm = 0                   otherwise
178      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
179      !!             tn  = ta
180      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
181      !!
182      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
183      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
184      !!----------------------------------------------------------------------
185      INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index
186      INTEGER , INTENT(in   )                               ::  kit000        ! first time-step index
187      INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers
188      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields
189      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields
190      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend
191      !!
192      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
193      REAL(wp) :: ztm, ztf     ! temporary scalars
194      !!----------------------------------------------------------------------
195
196      IF( kt == kit000 ) THEN
197         IF(lwp) WRITE(numout,*)
198         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
199         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
200      ENDIF
201      !
202      !                                              ! ----------------------- !
203      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
204         !                                           ! ----------------------- !
205         !
206         IF( neuler == 0 .AND. kt == kit000 ) THEN        ! Euler time-stepping at first time-step
207            !                                             ! (only swap)
208            DO jn = 1, kjpt
209               DO jk = 1, jpkm1                             
210                  DO jj = 1, jpj
211                     DO ji = 1, jpi
212                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptb <-- ptn
213                     END DO
214                  END DO
215               END DO
216            END DO
217         ELSE                                             ! general case (Leapfrog + Asselin filter
218            DO jn = 1, kjpt
219               DO jk = 1, jpkm1
220                  DO jj = 1, jpj
221                     DO ji = 1, jpi
222                        ztm = 0.25 * ( pta(ji,jj,jk,jn) + 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )  ! mean pt
223                        ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptn(ji,jj,jk,jn) )  ! Asselin filter on pt
224                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                      ! ptb <-- filtered ptn
225                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                            ! ptn <-- pta
226                        pta(ji,jj,jk,jn) = ztm                                                           ! pta <-- mean pt
227                     END DO
228                  END DO
229               END DO
230            END DO
231         ENDIF
232         !                                           ! ----------------------- !
233      ELSE                                           !    explicit hpg case    !
234         !                                           ! ----------------------- !
235         !
236         IF( neuler == 0 .AND. kt == kit000 ) THEN        ! Euler time-stepping at first time-step
237            DO jn = 1, kjpt
238               DO jk = 1, jpkm1
239                  DO jj = 1, jpj
240                     DO ji = 1, jpi
241                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                             ! ptn <-- pta
242                     END DO
243                  END DO
244               END DO
245            END DO
246         ELSE                                             ! general case (Leapfrog + Asselin filter
247            DO jn = 1, kjpt
248               DO jk = 1, jpkm1
249                  DO jj = 1, jpj
250                     DO ji = 1, jpi
251                        ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )       ! Asselin filter on t
252                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                     ! ptb <-- filtered ptn
253                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                           ! ptn <-- pta
254                     END DO
255                  END DO
256               END DO
257            END DO
258         ENDIF
259         !
260      ENDIF
261      !
262   END SUBROUTINE tra_nxt_fix
263
264   SUBROUTINE tra_nxt_vvl( kt, kit000,                    &
265      &                               ptb, ptn, pta, kjpt  )
266      !!----------------------------------------------------------------------
267      !!                   ***  ROUTINE tra_nxt_vvl  ***
268      !!
269      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
270      !!                and swap the tracer fields.
271      !!
272      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
273      !!              - save in (ta,sa) a thickness weighted average over the three
274      !!             time levels which will be used to compute rdn and thus the semi-
275      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
276      !!              - swap tracer fields to prepare the next time_step.
277      !!                This can be summurized for tempearture as:
278      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T
279      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )   
280      !!             ztm = 0                                otherwise
281      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
282      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
283      !!             tn  = ta
284      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
285      !!
286      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
287      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
288      !!----------------------------------------------------------------------
289      INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index
290      INTEGER , INTENT(in   )                               ::  kit000        ! first time-step index
291      INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers
292      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields
293      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields
294      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend
295      !!     
296      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices
297      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
298      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
299      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
300      !!----------------------------------------------------------------------
301
302      IF( kt == kit000 ) THEN
303         IF(lwp) WRITE(numout,*)
304         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
305         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
306      ENDIF
307      !
308      !                                              ! ----------------------- !
309      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
310         !                                           ! ----------------------- !
311         !
312         IF( neuler == 0 .AND. kt == kit000 ) THEN        ! Euler time-stepping at first time-step
313            DO jn = 1, kjpt                               ! (only swap)
314               DO jk = 1, jpkm1                             
315                  DO jj = 1, jpj
316                     DO ji = 1, jpi
317                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                    ! tn <-- ta
318                     END DO
319                  END DO
320               END DO
321            END DO
322         ELSE
323            DO jn = 1, kjpt                               ! (only swap)
324               DO jk = 1, jpkm1
325                  DO jj = 1, jpj
326                     DO ji = 1, jpi
327                        ze3t_b = fse3t_b(ji,jj,jk)
328                        ze3t_n = fse3t_n(ji,jj,jk)
329                        ze3t_a = fse3t_a(ji,jj,jk)
330                        !                                         ! tracer content at Before, now and after
331                        ztcb = ptb(ji,jj,jk,jn) *  ze3t_b 
332                        ztcn = ptn(ji,jj,jk,jn) *  ze3t_n 
333                        ztca = pta(ji,jj,jk,jn) *  ze3t_a 
334                        !
335                        !                                         ! Asselin filter on thickness and tracer content
336                        ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
337                        ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
338                        !
339                        !                                         ! filtered tracer including the correction
340                        ze3fr = 1.e0  / ( ze3t_n + ze3t_f )
341                        ztf   = ze3fr * ( ztcn   + ztc_f  )
342                        !                                         ! mean thickness and tracer
343                        ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b )
344                        ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   )
345                        !!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !!
346                        !!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr
347                        !                                         ! swap of arrays
348                        ptb(ji,jj,jk,jn) = ztf                             ! ptb <-- ptn + filter
349                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)              ! ptn <-- pta
350                        pta(ji,jj,jk,jn) = ztm                             ! pta <-- mean t
351                     END DO
352                  END DO
353               END DO
354            END DO
355         ENDIF
356         !                                           ! ----------------------- !
357      ELSE                                           !    explicit hpg case    !
358         !                                           ! ----------------------- !
359         !
360         IF( neuler == 0 .AND. kt == kit000 ) THEN        ! case of Euler time-stepping at first time-step
361            DO jn = 1, kjpt                               ! No filter nor thickness weighting computation required   
362               DO jk = 1, jpkm1                           ! ONLY swap                       
363                  DO jj = 1, jpj                             
364                     DO ji = 1, jpi
365                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! tn <-- ta
366                     END DO
367                  END DO
368               END DO
369            END DO
370            !                                             ! general case (Leapfrog + Asselin filter)
371         ELSE                                             ! apply filter on thickness weighted tracer and swap
372            DO jn = 1, kjpt     
373               DO jk = 1, jpkm1
374                  DO jj = 1, jpj
375                     DO ji = 1, jpi
376                        ze3t_b = fse3t_b(ji,jj,jk)
377                        ze3t_n = fse3t_n(ji,jj,jk)
378                        ze3t_a = fse3t_a(ji,jj,jk)
379                        !                                         ! tracer content at Before, now and after
380                        ztcb = ptb(ji,jj,jk,jn) *  ze3t_b 
381                        ztcn = ptn(ji,jj,jk,jn) *  ze3t_n   
382                        ztca = pta(ji,jj,jk,jn) *  ze3t_a 
383                        !
384                        !                                         ! Asselin filter on thickness and tracer content
385                        ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
386                        ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
387                        !
388                        !                                         ! filtered tracer including the correction
389                        ze3fr = 1.e0 / ( ze3t_n + ze3t_f )
390                        ztf   =  ( ztcn  + ztc_f ) * ze3fr
391                        !                                         ! swap of arrays
392                        ptb(ji,jj,jk,jn) = ztf                  ! tb <-- tn filtered
393                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)   ! tn <-- ta
394                     END DO
395                  END DO
396               END DO
397            END DO
398         ENDIF
399      ENDIF
400      !
401   END SUBROUTINE tra_nxt_vvl
402
403   !!======================================================================
404END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.