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

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/tranxt.F90 @ 2120

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

apply cosmetic changes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 19.2 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 obc_oce
35   USE obctra          ! open boundary condition (obc_tra routine)
36   USE bdytra          ! Unstructured open boundary condition (bdy_tra routine)
37   USE in_out_manager  ! I/O manager
38   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
39   USE prtctl          ! Print control
40   USE traswp          ! swap array
41   USE agrif_opa_update
42   USE agrif_opa_interp
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
109#if defined key_obc
110      IF( lk_obc )   CALL obc_tra( kt )  ! OBC open boundaries
111#endif
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#if defined key_obc || defined key_bdy || defined key_agrif
120      CALL tra_swap
121#endif
122 
123      ! set time step size (Euler/Leapfrog)
124      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt(:) =     rdttra(:)      ! at nit000             (Euler)
125      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
126      ENDIF
127
128      ! trends computation initialisation
129      IF( l_trdtra )   THEN                    !* store now fields before applying the Asselin filter
130         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
131         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsn(:,:,:,jp_sal)
132      ENDIF
133
134      ! Leap-Frog + Asselin filter time stepping
135      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
136      ELSE                  ;   CALL tra_nxt_fix( kt, tsb, tsn, tsa, jpts )  ! fixed    volume level
137      ENDIF
138
139#if defined key_agrif
140      CALL tra_unswap
141      ! Update tracer at AGRIF zoom boundaries
142      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
143      CALL tra_swap
144#endif     
145
146      ! trends computation
147      IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt     
148         DO jk = 1, jpkm1
149            zfact = 1.e0 / r2dt(jk)             
150            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
151            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
152         END DO
153         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt )
154         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds )
155         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
156      END IF
157
158      !                        ! control print
159      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
160         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
161      !
162   END SUBROUTINE tra_nxt
163
164
165   SUBROUTINE tra_nxt_fix( kt, ptb, ptn, pta, kjpt )
166      !!----------------------------------------------------------------------
167      !!                   ***  ROUTINE tra_nxt_fix  ***
168      !!
169      !! ** Purpose :   fixed volume: apply the Asselin time filter and
170      !!                swap the tracer fields.
171      !!
172      !! ** Method  : - Apply a Asselin time filter on now fields.
173      !!              - save in (ta,sa) an average over the three time levels
174      !!             which will be used to compute rdn and thus the semi-implicit
175      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
176      !!              - swap tracer fields to prepare the next time_step.
177      !!                This can be summurized for tempearture as:
178      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T
179      !!             ztm = 0                   otherwise
180      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
181      !!             tn  = ta
182      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
183      !!
184      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
185      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
186      !!----------------------------------------------------------------------
187      INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index
188      INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers
189      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields
190      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields
191      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend
192      !!
193      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
194      REAL(wp) :: ztm, ztf     ! temporary scalars
195      !!----------------------------------------------------------------------
196
197      IF( kt == nit000 )  THEN
198         IF(lwp) WRITE(numout,*)
199         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
200         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
201      ENDIF
202      !
203      !
204      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
205         !                                             ! (only swap)
206         DO jn = 1, kjpt
207            DO jk = 1, jpkm1                             
208               ptn(:,:,jk,jn) = pta(:,:,jk,jn)     ! ptb <-- ptn
209            END DO
210         END DO
211         !                                             
212      ELSE                                           ! general case (Leapfrog + Asselin filter
213         !
214         !                                           ! ----------------------- !
215         IF( ln_dynhpg_imp ) THEN                    ! semi-implicite hpg case !
216            !                                        ! ----------------------- !
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            !                                           ! ----------------------- !
231         ELSE                                           !    explicit hpg case    !
232            !                                           ! ----------------------- !
233            DO jn = 1, kjpt
234               DO jk = 1, jpkm1
235                  DO jj = 1, jpj
236                     DO ji = 1, jpi
237                        ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )       ! Asselin filter on t
238                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                     ! ptb <-- filtered ptn
239                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                           ! ptn <-- pta
240                     END DO
241                  END DO
242               END DO
243            END DO
244            !
245         ENDIF
246         !
247      ENDIF
248      !
249   END SUBROUTINE tra_nxt_fix
250
251
252   SUBROUTINE tra_nxt_vvl( kt, ptb, ptn, pta, kjpt )
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      INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers
278      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields
279      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields
280      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend
281      !!     
282      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices
283      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
284      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
285      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
286      !!----------------------------------------------------------------------
287
288      IF( kt == nit000 ) THEN
289         IF(lwp) WRITE(numout,*)
290         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
291         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
292      ENDIF
293      !
294      !
295      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
296         !                                             ! (only swap)
297         DO jn = 1, kjpt
298            DO jk = 1, jpkm1                             
299               ptn(:,:,jk,jn) = pta(:,:,jk,jn)         ! ptb <-- ptn
300            END DO
301         END DO
302         !                                             
303      ELSE                                              ! general case (Leapfrog + Asselin filter
304         !
305         !                                              ! ----------------------- !
306         IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
307            !                                           ! ----------------------- !
308            DO jn = 1, kjpt                               
309               DO jk = 1, jpkm1
310                  DO jj = 1, jpj
311                     DO ji = 1, jpi
312                        ze3t_b = fse3t_b(ji,jj,jk)
313                        ze3t_n = fse3t_n(ji,jj,jk)
314                        ze3t_a = fse3t_a(ji,jj,jk)
315                        !                                         ! tracer content at Before, now and after
316                        ztcb = ptb(ji,jj,jk,jn) *  ze3t_b 
317                        ztcn = ptn(ji,jj,jk,jn) *  ze3t_n 
318                        ztca = pta(ji,jj,jk,jn) *  ze3t_a 
319                        !
320                        !                                         ! Asselin filter on thickness and tracer content
321                        ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
322                        ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
323                        !
324                        !                                         ! filtered tracer including the correction
325                        ze3fr = 1.e0  / ( ze3t_n + ze3t_f )
326                        ztf   = ze3fr * ( ztcn   + ztc_f  )
327                        !                                         ! mean thickness and tracer
328                        ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b )
329                        ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   )
330                        !!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !!
331                        !!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr
332                        !                                         ! swap of arrays
333                        ptb(ji,jj,jk,jn) = ztf                             ! ptb <-- ptn + filter
334                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)              ! ptn <-- pta
335                        pta(ji,jj,jk,jn) = ztm                             ! pta <-- mean t
336                     END DO
337                  END DO
338               END DO
339            END DO
340            !                                        ! ----------------------- !
341         ELSE                                        !    explicit hpg case    !
342            !                                        ! ----------------------- !
343            DO jn = 1, kjpt     
344               DO jk = 1, jpkm1
345                  DO jj = 1, jpj
346                     DO ji = 1, jpi
347                        ze3t_b = fse3t_b(ji,jj,jk)
348                        ze3t_n = fse3t_n(ji,jj,jk)
349                        ze3t_a = fse3t_a(ji,jj,jk)
350                        !                                         ! tracer content at Before, now and after
351                        ztcb = ptb(ji,jj,jk,jn) *  ze3t_b 
352                        ztcn = ptn(ji,jj,jk,jn) *  ze3t_n   
353                        ztca = pta(ji,jj,jk,jn) *  ze3t_a 
354                        !
355                        !                                         ! Asselin filter on thickness and tracer content
356                        ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
357                        ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
358                        !
359                        !                                         ! filtered tracer including the correction
360                        ze3fr = 1.e0 / ( ze3t_n + ze3t_f )
361                        ztf   =  ( ztcn  + ztc_f ) * ze3fr
362                        !                                         ! swap of arrays
363                        ptb(ji,jj,jk,jn) = ztf                    ! tb <-- tn filtered
364                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)       ! tn <-- ta
365                     END DO
366                  END DO
367               END DO
368            END DO
369            !
370         ENDIF
371         !
372      ENDIF
373      !
374   END SUBROUTINE tra_nxt_vvl
375
376   !!======================================================================
377END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.