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

source: branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90 @ 2068

Last change on this file since 2068 was 2068, checked in by mlelod, 14 years ago

ticket: #663 ensuring restartability and conservation
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 17.4 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-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
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 sbc_oce         ! surface boundary condition: ocean
28   USE zdf_oce         ! ???
29   USE domvvl          ! variable volume
30   USE dynspg_oce      ! surface     pressure gradient variables
31   USE dynhpg          ! hydrostatic pressure gradient
32   USE trdmod_oce      ! ocean variables trends
33   USE trdmod          ! ocean active tracers trends
34   USE phycst
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 traqsr          ! penetrative solar radiation (needed for nksr)
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
49   REAL(wp)                 ::   rbcp            ! Brown & Campana parameters for semi-implicit hpg
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.3 , LOCEAN-IPSL (2010)
56   !! $Id$
57   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_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      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
87      !!----------------------------------------------------------------------
88      USE oce, ONLY :    ztrdt => ua   ! use ua as 3D workspace   
89      USE oce, ONLY :    ztrds => va   ! use va as 3D workspace   
90      !!
91      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
92      !!
93      INTEGER  ::   jk      ! dummy loop indices
94      REAL(wp) ::   zfact   ! temporary scalars
95      !!----------------------------------------------------------------------
96
97      IF( kt == nit000 ) THEN          !==  initialisation  ==!
98         IF(lwp) WRITE(numout,*)
99         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
100         IF(lwp) WRITE(numout,*) '~~~~~~~'
101         !
102         rbcp    = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)       ! Brown & Campana parameter for semi-implicit hpg
103      ENDIF
104      !                                      ! set time step size (Euler/Leapfrog)
105      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
106      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
107      ENDIF
108
109      !                                !==  Update after tracer on domain lateral boundaries  ==!
110      !
111      CALL lbc_lnk( ta, 'T', 1. )            ! local domain boundaries  (T-point, unchanged sign)
112      CALL lbc_lnk( sa, 'T', 1. )
113      !
114#if defined key_obc
115      CALL obc_tra( kt )                     ! OBC open boundaries
116#endif
117#if defined key_bdy
118      CALL bdy_tra( kt )                     ! BDY open boundaries
119#endif
120#if defined key_agrif
121      CALL Agrif_tra                         ! AGRIF zoom boundaries
122#endif
123 
124      IF( l_trdtra ) THEN              ! trends computation: store now fields before applying the Asselin filter
125         ztrdt(:,:,:) = tn(:,:,:)     
126         ztrds(:,:,:) = sn(:,:,:)
127      ENDIF
128
129      !                                !==  modifed Leap-Frog + Asselin filter (modified LF-RA) scheme  ==!
130      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt )      ! variable volume level (vvl)
131      ELSE                  ;   CALL tra_nxt_fix( kt )      ! fixed    volume level
132      ENDIF
133
134#if defined key_agrif
135      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )   ! Update tracer at AGRIF zoom boundaries (children only)
136#endif     
137
138      IF( l_trdtra ) THEN              ! trends computation: 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 temperature as:
167      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
168      !!             ztm = 0                                   otherwise
169      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
170      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
171      !!             tn  = ta
172      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
173      !!
174      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
175      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
176      !!----------------------------------------------------------------------
177      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
178      !!
179      INTEGER  ::   ji, jj, jk      ! dummy loop indices
180      REAL(wp) ::   zt_d, zs_d      ! temporary scalars
181      REAL(wp) ::   ztn, zsn        !    -         -
182      !!----------------------------------------------------------------------
183
184      IF( kt == nit000 ) THEN
185         IF(lwp) WRITE(numout,*)
186         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
187         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
188      ENDIF
189      !
190      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
191         DO jk = 1, jpkm1                                                 ! (only swap)
192            tn(:,:,jk) = ta(:,:,jk)                                       ! tn <-- ta
193            sn(:,:,jk) = sa(:,:,jk)                                       ! sn <-- sa
194         END DO
195      ELSE                                             ! General case (Leapfrog + Asselin filter)
196         DO jk = 1, jpkm1
197            DO jj = 1, jpj
198               DO ji = 1, jpi
199                  IF( ln_dynhpg_imp ) THEN                  ! implicit hpg: keep tn, sn in memory
200                     ztn = tn(ji,jj,jk)
201                     zsn = sn(ji,jj,jk)
202                  ENDIF
203                  !                                         ! time laplacian on tracers
204                  !                                         ! used for both Asselin and Brown & Campana filters
205                  zt_d = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk)
206                  zs_d = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk)
207                  !
208                  !                                         ! swap of arrays
209                  tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_d               ! tb <-- tn filtered
210                  sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_d               ! sb <-- sn filtered
211                  tn(ji,jj,jk) = ta(ji,jj,jk)                             ! tn <-- ta
212                  sn(ji,jj,jk) = sa(ji,jj,jk)                             ! sn <-- sa
213                  !                                         ! semi imlicit hpg computation (Brown & Campana)
214                  IF( ln_dynhpg_imp ) THEN
215                     ta(ji,jj,jk) = ztn + rbcp * zt_d                     ! ta <-- Brown & Campana average
216                     sa(ji,jj,jk) = zsn + rbcp * zs_d                     ! sa <-- Brown & Campana average
217                  ENDIF
218               END DO
219            END DO
220         END DO
221      ENDIF
222      !
223   END SUBROUTINE tra_nxt_fix
224
225
226   SUBROUTINE tra_nxt_vvl( kt )
227      !!----------------------------------------------------------------------
228      !!                   ***  ROUTINE tra_nxt_vvl  ***
229      !!
230      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
231      !!                and swap the tracer fields.
232      !!
233      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
234      !!              - save in (ta,sa) a thickness weighted average over the three
235      !!             time levels which will be used to compute rdn and thus the semi-
236      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
237      !!              - swap tracer fields to prepare the next time_step.
238      !!                This can be summurized for temperature as:
239      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
240      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
241      !!             ztm = 0                                                       otherwise
242      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
243      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
244      !!             tn  = ta
245      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
246      !!
247      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
248      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
249      !!----------------------------------------------------------------------
250      INTEGER, INTENT(in) ::   kt                  ! ocean time-step index
251      !!     
252      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
253      REAL     ::   ze3t_a, ze3t_n, ze3t_b         ! temporary scalars
254      REAL     ::   ztc_a, ztc_n, ztc_b            !    -         -
255      REAL     ::   zsc_a, zsc_n, zsc_b            !    -         -
256      REAL     ::   ztc_f, zsc_f, ztc_d, zsc_d     !    -         -
257      REAL     ::   ze3t_f, ze3t_d                 !    -         -
258      REAL     ::   zfact1, zfact2                 !    -         -
259      !!----------------------------------------------------------------------
260     
261      IF( kt == nit000 ) THEN
262         IF(lwp) WRITE(numout,*)
263         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
264         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
265      ENDIF
266
267      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
268         !                                                  ! (only swap)
269         DO jk = 1, jpkm1                                   ! tn <-- ta
270            tn(:,:,jk) = ta(:,:,jk)                         ! sn <-- sa
271            sn(:,:,jk) = sa(:,:,jk)
272         END DO
273         !                                             ! General case (Leapfrog + Asselin filter)
274      ELSE                                             ! apply filter on thickness weighted tracer and swap
275         DO jk = 1, jpkm1
276            zfact1 = atfp * rdttra(jk)
277            zfact2 = zfact1 / rau0
278            DO jj = 1, jpj
279               DO ji = 1, jpi
280                  !                                         ! scale factors at Before, now and after
281                  ze3t_b = fse3t_b(ji,jj,jk)
282                  ze3t_n = fse3t_n(ji,jj,jk)
283                  ze3t_a = fse3t_a(ji,jj,jk)
284                  ze3t_d = fse3t_d(ji,jj,jk)
285                  !                                         ! tracer content at Before, now and after
286                  ztc_b  = tb(ji,jj,jk) * ze3t_b   ;   zsc_b = sb(ji,jj,jk) * ze3t_b
287                  ztc_n  = tn(ji,jj,jk) * ze3t_n   ;   zsc_n = sn(ji,jj,jk) * ze3t_n
288                  ztc_a  = ta(ji,jj,jk) * ze3t_a   ;   zsc_a = sa(ji,jj,jk) * ze3t_a
289                  !
290                  !                                         ! Time laplacian on tracer contents
291                  !                                         ! used for both Asselin and Brown & Campana filters
292                  ztc_d  = ztc_a - 2. * ztc_n + ztc_b
293                  zsc_d  = zsc_a - 2. * zsc_n + zsc_b
294                  !                                         ! Asselin Filter on thicknesses and tracer contents
295                  ze3t_f = ze3t_n + atfp * ze3t_d
296                  ztc_f  = ztc_n  + atfp * ztc_d
297                  zsc_f  = zsc_n  + atfp * zsc_d
298                  !                                         ! Filter correction
299                  IF( jk == 1 ) THEN
300                     ! WRITE(numout,*) 'filter correction: sbc_trd_hc_n'
301                     ze3t_f = ze3t_f - zfact2 * ( emp_b       (ji,jj) - emp         (ji,jj) )
302                     ztc_f  = ztc_f  - zfact1 * ( sbc_trd_hc_n(ji,jj) - sbc_trd_hc_b(ji,jj) )
303                  ENDIF
304                  IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN
305                     ! WRITE(numout,*) 'jk =', jk
306                     ! WRITE(numout,*) 'filter correction: qsr_trd_hc_n'
307                     ztc_f  = ztc_f  - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) )
308                  ENDIF
309                                                          ! swap of arrays
310                  ze3t_f = 1.e0 / ze3t_f
311                  tb(ji,jj,jk) = ztc_f * ze3t_f                           ! tb <-- tn filtered
312                  sb(ji,jj,jk) = zsc_f * ze3t_f                           ! sb <-- sn filtered
313                  tn(ji,jj,jk) = ta(ji,jj,jk)                             ! tn <-- ta
314                  sn(ji,jj,jk) = sa(ji,jj,jk)                             ! sn <-- sa
315                  !                                         ! semi imlicit hpg computation (Brown & Campana)
316                  IF( ln_dynhpg_imp ) THEN
317                     ze3t_d       = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
318                     ta(ji,jj,jk) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
319                     sa(ji,jj,jk) = ze3t_d * ( zsc_n  + rbcp * zsc_d  )   ! sa <-- Brown & Campana average
320                  ENDIF
321               END DO
322            END DO
323         END DO
324      ENDIF
325      !
326   END SUBROUTINE tra_nxt_vvl
327
328   !!======================================================================
329END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.