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

source: branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90 @ 1361

Last change on this file since 1361 was 1361, checked in by rblod, 15 years ago

dev004_VVL: new organisation see ticket #389

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.0 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   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   tra_nxt       : time stepping on temperature and salinity
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers variables
22   USE dom_oce         ! ocean space and time domain variables
23   USE zdf_oce         ! ???
24   USE dynspg_oce      ! surface pressure gradient variables
25   USE trdmod_oce      ! ocean variables trends
26   USE trdmod          ! ocean active tracers trends
27   USE phycst
28   USE obctra          ! open boundary condition (obc_tra routine)
29   USE bdytra          ! Unstructured open boundary condition (bdy_tra routine)
30   USE in_out_manager  ! I/O manager
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE prtctl          ! Print control
33   USE agrif_opa_update
34   USE agrif_opa_interp
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   tra_nxt    ! routine called by step.F90
40
41   REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
47   !! $Id$
48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE tra_nxt( kt )
54      !!----------------------------------------------------------------------
55      !!                   ***  ROUTINE tranxt  ***
56      !!
57      !! ** Purpose :   Apply the boundary condition on the after temperature 
58      !!             and salinity fields, achieved the time stepping by adding
59      !!             the Asselin filter on now fields and swapping the fields.
60      !!
61      !! ** Method  :   At this stage of the computation, ta and sa are the
62      !!             after temperature and salinity as the time stepping has
63      !!             been performed in trazdf_imp or trazdf_exp module.
64      !!
65      !!              - Apply lateral boundary conditions on (ta,sa)
66      !!             at the local domain   boundaries through lbc_lnk call,
67      !!             at the radiative open boundaries (lk_obc=T),
68      !!             at the relaxed   open boundaries (lk_bdy=T), and
69      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
70      !!
71      !!              - Update lateral boundary conditions on AGRIF children
72      !!             domains (lk_agrif=T)
73      !!
74      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
75      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
76      !!----------------------------------------------------------------------
77      USE oce, ONLY :    ztrdt => ua   ! use ua as 3D workspace   
78      USE oce, ONLY :    ztrds => va   ! use va as 3D workspace   
79      !!
80      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
81      !!
82      INTEGER  ::   jk    ! dummy loop indices
83      REAL(wp) ::   zfact ! temporary scalars
84      !!----------------------------------------------------------------------
85
86      IF( kt == nit000 ) THEN
87         IF(lwp) WRITE(numout,*)
88         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
89         IF(lwp) WRITE(numout,*) '~~~~~~~'
90      ENDIF
91
92      ! Update after tracer on domain lateral boundaries
93      !
94      CALL lbc_lnk( ta, 'T', 1. )      ! local domain boundaries  (T-point, unchanged sign)
95      CALL lbc_lnk( sa, 'T', 1. )
96      !
97#if defined key_obc
98      CALL obc_tra( kt )               ! OBC open boundaries
99#endif
100#if defined key_bdy
101      CALL bdy_tra( kt )               ! BDY open boundaries
102#endif
103#if defined key_agrif
104      CALL Agrif_tra                   ! AGRIF zoom boundaries
105#endif
106 
107      ! set time step size (Euler/Leapfrog)
108      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
109      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
110      ENDIF
111
112      ! trends computation initialisation
113      IF( l_trdtra ) THEN              ! store now fields before applying the Asselin filter
114         ztrdt(:,:,:) = tn(:,:,:)     
115         ztrds(:,:,:) = sn(:,:,:)
116      ENDIF
117
118      ! Leap-Frog + Asselin filter time stepping
119      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt )      ! variable volume level (vvl)
120      ELSE                  ;   CALL tra_nxt_fix( kt )      ! fixed    volume level
121      ENDIF
122
123#if defined key_agrif
124      ! Update tracer at AGRIF zoom boundaries
125      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
126#endif     
127
128      ! trends diagnostics :  Asselin filter trend : (tb filtered - tb)/2dt
129      IF( l_trdtra ) THEN     
130         DO jk = 1, jpkm1
131            zfact = 1.e0 / r2dt_t(jk)             
132            ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact
133            ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact
134         END DO
135         CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt )
136      END IF
137      !                        ! control print
138      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   &
139         &                       tab3d_2=sn, clinfo2=       ' Sn: ', mask2=tmask )
140      !
141   END SUBROUTINE tra_nxt
142
143
144   SUBROUTINE tra_nxt_fix( kt )
145      !!----------------------------------------------------------------------
146      !!                   ***  ROUTINE tra_nxt_fix  ***
147      !!
148      !! ** Purpose :   fixed volume: apply the Asselin time filter and
149      !!                swap the tracer fields.
150      !!
151      !! ** Method  : - Apply a Asselin time filter on now fields.
152      !!              - save in (ta,sa) an average over the three time levels
153      !!             which will be used to compute rdn and thus the semi-implicit
154      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
155      !!              - swap tracer fields to prepare the next time_step.
156      !!                This can be summurized for tempearture as:
157      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T
158      !!             ztm = 0                   otherwise
159      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
160      !!             tn  = ta
161      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
162      !!
163      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
164      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
165      !!----------------------------------------------------------------------
166      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
167      !!
168      INTEGER  ::   ji, jj, jk   ! dummy loop indices
169      REAL(wp) ::   ztm, ztf     ! temporary scalars
170      REAL(wp) ::   zsm, zsf     !    -         -
171      !!----------------------------------------------------------------------
172
173      IF( kt == nit000 ) THEN
174         IF(lwp) WRITE(numout,*)
175         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
176         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
177      ENDIF
178      !
179      !                                              ! ----------------------- !
180      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
181         !                                           ! ----------------------- !
182         !
183         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
184            DO jk = 1, jpkm1
185               DO jj = 1, jpj
186                  DO ji = 1, jpi
187                     ztm = 0.25 * ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) )      ! mean t
188                     zsm = 0.25 * ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) )
189                     tb(ji,jj,jk) = tn(ji,jj,jk)                                           ! tb <-- tn
190                     sb(ji,jj,jk) = sn(ji,jj,jk)
191                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tb <-- tn
192                     sn(ji,jj,jk) = sa(ji,jj,jk)
193                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t
194                     sa(ji,jj,jk) = zsm
195                  END DO
196               END DO
197            END DO
198         ELSE                                             ! general case (Leapfrog + Asselin filter
199            DO jk = 1, jpkm1
200               DO jj = 1, jpj
201                  DO ji = 1, jpi
202                     ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! mean t
203                     zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
204                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t
205                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
206                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn
207                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
208                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
209                     sn(ji,jj,jk) = sa(ji,jj,jk)
210                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t
211                     sa(ji,jj,jk) = zsm
212                  END DO
213               END DO
214            END DO
215         ENDIF
216         !                                           ! ----------------------- !
217      ELSE                                           !    explicit hpg case    !
218         !                                           ! ----------------------- !
219         !
220         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
221            DO jk = 1, jpkm1
222               DO jj = 1, jpj
223                  DO ji = 1, jpi
224                     tb(ji,jj,jk) = tn(ji,jj,jk)                                           ! tb <-- tn
225                     sb(ji,jj,jk) = sn(ji,jj,jk)
226                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
227                     sn(ji,jj,jk) = sa(ji,jj,jk)
228                  END DO
229               END DO
230            END DO
231         ELSE                                             ! general case (Leapfrog + Asselin filter
232            DO jk = 1, jpkm1
233               DO jj = 1, jpj
234                  DO ji = 1, jpi
235!RBvvl for reproducibility
236!                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t
237!                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
238!                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn
239!                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
240                     tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
241                     sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
242                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
243                     sn(ji,jj,jk) = sa(ji,jj,jk)
244                  END DO
245               END DO
246            END DO
247         ENDIF
248      ENDIF
249      !
250   END SUBROUTINE tra_nxt_fix
251
252
253   SUBROUTINE tra_nxt_vvl( kt )
254      !!----------------------------------------------------------------------
255      !!                   ***  ROUTINE tra_nxt_vvl  ***
256      !!
257      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
258      !!                and swap the tracer fields.
259      !!
260      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
261      !!              - save in (ta,sa) a thickness weighted average over the three
262      !!             time levels which will be used to compute rdn and thus the semi-
263      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
264      !!              - swap tracer fields to prepare the next time_step.
265      !!                This can be summurized for tempearture as:
266      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T
267      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )   
268      !!             ztm = 0                                otherwise
269      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
270      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
271      !!             tn  = ta
272      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
273      !!
274      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
275      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
276      !!----------------------------------------------------------------------
277      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
278      !!     
279     
280      ! Not yet ready
281      WRITE(*,*) 'tra_next_vvl : you should not be there'
282      STOP
283      !
284   END SUBROUTINE tra_nxt_vvl
285
286   !!======================================================================
287END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.