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

source: branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/TRA/tranxt.F90 @ 2208

Last change on this file since 2208 was 2208, checked in by rblod, 14 years ago

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