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

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

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