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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA/tranxt.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

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