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 @ 1975

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

ticket: #663 MLF: first part
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 17.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-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_m, zs_m      ! 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_m = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk)
206                  zs_m = 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_m               ! tb <-- tn filtered
210                  sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_m               ! 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_m                     ! ta <-- Brown & Campana average
216                     sa(ji,jj,jk) = zsn + rbcp * zs_m                     ! 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_m, zsc_m     !    -         -
257      REAL     ::   ze3t_f, ze3t_m                 !    -         -
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 * r2dt_t(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                  !                                         ! tracer content at Before, now and after
285                  ztc_b  = tb(ji,jj,jk) * ze3t_b   ;   zsc_b = sb(ji,jj,jk) * ze3t_b
286                  ztc_n  = tn(ji,jj,jk) * ze3t_n   ;   zsc_n = sn(ji,jj,jk) * ze3t_n
287                  ztc_a  = ta(ji,jj,jk) * ze3t_a   ;   zsc_a = sa(ji,jj,jk) * ze3t_a
288                  !
289                  !                                         ! Time laplacian on thickness and tracer content
290                  !                                         ! used for both Asselin and Brown & Campana filters
291                  ze3t_m = ze3t_a - 2. * ze3t_n + ze3t_b
292                  ztc_m  = ztc_a  - 2. * ztc_n  + ztc_b
293                  zsc_m  = zsc_a  - 2. * zsc_n  + zsc_b
294                  !                                         ! Asselin Filter + correction
295                  ze3t_f = ze3t_n + atfp * ze3t_m
296                  ztc_f  = ztc_n  + atfp * ztc_m
297                  zsc_f  = zsc_n  + atfp * zsc_m
298
299                  IF( jk == 1 ) THEN
300                     ze3t_f = ze3t_f - zfact2 * ( emp_b       (ji,jj)    - emp         (ji,jj)    )
301                     ztc_f  = ztc_f  - zfact1 * ( sbc_trd_hc_n(ji,jj)    - sbc_trd_hc_b(ji,jj)    )
302                  ENDIF
303                  IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN
304                     ztc_f  = ztc_f  - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) )
305                  ENDIF
306                  !                                         ! swap of arrays
307                  ze3t_f = 1.e0 / ze3t_f
308                  tb(ji,jj,jk) = ztc_f * ze3t_f                           ! tb <-- tn filtered
309                  sb(ji,jj,jk) = zsc_f * ze3t_f                           ! sb <-- sn filtered
310                  tn(ji,jj,jk) = ta(ji,jj,jk)                             ! tn <-- ta
311                  sn(ji,jj,jk) = sa(ji,jj,jk)                             ! sn <-- sa
312                  !                                         ! semi imlicit hpg computation (Brown & Campana)
313                  IF( ln_dynhpg_imp ) THEN
314                     ze3t_m = 1.e0 / ( ze3t_n + rbcp * ze3t_m )
315                     ta(ji,jj,jk) =  ( ztc_n  + rbcp * ztc_m  ) * ze3t_m  ! ta <-- Brown & Campana average
316                     sa(ji,jj,jk) =  ( zsc_n  + rbcp * zsc_m  ) * ze3t_m  ! sa <-- Brown & Campana average
317                  ENDIF
318               END DO
319            END DO
320         END DO
321      ENDIF
322      !
323   END SUBROUTINE tra_nxt_vvl
324
325   !!======================================================================
326END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.