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

source: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 2007

Last change on this file since 2007 was 2007, checked in by smasson, 14 years ago

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