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

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

merge LOCEAN 2010 developments branches
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 21.5 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   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   tra_nxt       : time stepping on tracers
23   !!   tra_nxt_fix   : time stepping on tracers : fixed    volume case
24   !!   tra_nxt_vvl   : time stepping on tracers : variable volume case
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers variables
27   USE dom_oce         ! ocean space and time domain variables
28   USE sbc_oce         ! surface boundary condition: ocean
29   USE zdf_oce         ! ???
30   USE domvvl          ! variable volume
31   USE dynspg_oce      ! surface     pressure gradient variables
32   USE dynhpg          ! hydrostatic pressure gradient
33   USE trdmod_oce      ! ocean space and time domain variables
34   USE trdtra          ! ocean active tracers trends
35   USE phycst
36   USE obc_oce
37   USE obctra          ! open boundary condition (obc_tra routine)
38   USE bdytra          ! Unstructured open boundary condition (bdy_tra routine)
39   USE in_out_manager  ! I/O manager
40   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
41   USE prtctl          ! Print control
42   USE traqsr          ! penetrative solar radiation (needed for nksr)
43   USE traswp          ! swap array
44   USE agrif_opa_update
45   USE agrif_opa_interp
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   tra_nxt       ! routine called by step.F90
51   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
52   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
53
54   REAL(wp)                 ::   rbcp            ! Brown & Campana parameters for semi-implicit hpg
55   REAL(wp), DIMENSION(jpk) ::   r2dt   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
56
57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64
65CONTAINS
66
67   SUBROUTINE tra_nxt( kt )
68      !!----------------------------------------------------------------------
69      !!                   ***  ROUTINE tranxt  ***
70      !!
71      !! ** Purpose :   Apply the boundary condition on the after temperature 
72      !!             and salinity fields, achieved the time stepping by adding
73      !!             the Asselin filter on now fields and swapping the fields.
74      !!
75      !! ** Method  :   At this stage of the computation, ta and sa are the
76      !!             after temperature and salinity as the time stepping has
77      !!             been performed in trazdf_imp or trazdf_exp module.
78      !!
79      !!              - Apply lateral boundary conditions on (ta,sa)
80      !!             at the local domain   boundaries through lbc_lnk call,
81      !!             at the radiative open boundaries (lk_obc=T),
82      !!             at the relaxed   open boundaries (lk_bdy=T), and
83      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
84      !!
85      !!              - Update lateral boundary conditions on AGRIF children
86      !!             domains (lk_agrif=T)
87      !!
88      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
89      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
90      !!----------------------------------------------------------------------
91      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
92      !!
93      INTEGER  ::   jk    ! dummy loop indices
94      REAL(wp) ::   zfact ! local scalars
95      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
96      !!----------------------------------------------------------------------
97
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
101         IF(lwp) WRITE(numout,*) '~~~~~~~'
102         !
103         rbcp    = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)       ! Brown & Campana parameter for semi-implicit hpg
104      ENDIF
105
106      ! Update after tracer on domain lateral boundaries
107      !
108      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )      ! local domain boundaries  (T-point, unchanged sign)
109      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
110      !
111#if defined key_obc || defined key_bdy || defined key_agrif
112      CALL tra_unswap
113#endif
114
115#if defined key_obc 
116      IF( lk_obc )   CALL obc_tra( kt )  ! OBC open boundaries
117#endif
118#if defined key_bdy 
119      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
120#endif
121#if defined key_agrif
122      CALL Agrif_tra                   ! AGRIF zoom boundaries
123#endif
124
125#if defined key_obc || defined key_bdy || defined key_agrif
126      CALL tra_swap
127#endif
128 
129      ! set time step size (Euler/Leapfrog)
130      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt(:) =     rdttra(:)      ! at nit000             (Euler)
131      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
132      ENDIF
133
134      ! trends computation initialisation
135      IF( l_trdtra )   THEN                    !* store now fields before applying the Asselin filter
136         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
137         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsn(:,:,:,jp_sal)
138      ENDIF
139
140      ! Leap-Frog + Asselin filter time stepping
141      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
142      ELSE                  ;   CALL tra_nxt_fix( kt,        tsb, tsn, tsa, jpts )  ! fixed    volume level
143      ENDIF
144
145#if defined key_agrif
146      CALL tra_unswap
147      ! Update tracer at AGRIF zoom boundaries
148      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
149      CALL tra_swap
150#endif     
151
152      ! trends computation
153      IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt     
154         DO jk = 1, jpkm1
155            zfact = 1.e0 / r2dt(jk)             
156            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
157            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
158         END DO
159         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt )
160         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds )
161         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
162      END IF
163
164      !                        ! control print
165      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
166         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
167      !
168   END SUBROUTINE tra_nxt
169
170
171   SUBROUTINE tra_nxt_fix( kt, ptb, ptn, pta, kjpt )
172      !!----------------------------------------------------------------------
173      !!                   ***  ROUTINE tra_nxt_fix  ***
174      !!
175      !! ** Purpose :   fixed volume: apply the Asselin time filter and
176      !!                swap the tracer fields.
177      !!
178      !! ** Method  : - Apply a Asselin time filter on now fields.
179      !!              - save in (ta,sa) an average over the three time levels
180      !!             which will be used to compute rdn and thus the semi-implicit
181      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
182      !!              - swap tracer fields to prepare the next time_step.
183      !!                This can be summurized for tempearture as:
184      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
185      !!             ztm = 0                                   otherwise
186      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
187      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
188      !!             tn  = ta 
189      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
190      !!
191      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
192      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
193      !!----------------------------------------------------------------------
194      INTEGER , INTENT(in   )                               ::  kt       ! ocean time-step index
195      INTEGER , INTENT(in   )                               ::  kjpt     ! number of tracers
196      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
197      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
198      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
199      !!
200      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
201      REAL(wp) :: ztd, ztm         ! temporary scalars
202      !!----------------------------------------------------------------------
203
204      IF( kt == nit000 )  THEN
205         IF(lwp) WRITE(numout,*)
206         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
207         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
208      ENDIF
209      !
210      !
211      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
212         !                                             ! (only swap)
213         DO jn = 1, kjpt
214            DO jk = 1, jpkm1
215               ptn(:,:,jk,jn) = pta(:,:,jk,jn)     ! ptb <-- ptn
216            END DO
217         END DO
218         !                                             
219      ELSE                                           ! general case (Leapfrog + Asselin filter
220         !
221         !                                           ! ----------------------- !
222         IF( ln_dynhpg_imp ) THEN                    ! semi-implicite hpg case !
223            !                                        ! ----------------------- !
224            DO jn = 1, kjpt
225               DO jk = 1, jpkm1
226                  DO jj = 1, jpj
227                     DO ji = 1, jpi
228                        ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    !  time laplacian on tracers
229                        ztm = ptn(ji,jj,jk,jn) + rbcp * ztd                                 !  used for both Asselin and Brown & Campana filters
230                        !
231                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn
232                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta
233                        pta(ji,jj,jk,jn) = ztm                                              ! pta <-- Brown & Campana average
234                     END DO
235                  END DO
236               END DO
237            END DO
238            !                                           ! ----------------------- !
239         ELSE                                           !    explicit hpg case    !
240            !                                           ! ----------------------- !
241            DO jn = 1, kjpt
242               DO jk = 1, jpkm1
243                  DO jj = 1, jpj
244                     DO ji = 1, jpi
245                        ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    ! time laplacian on tracers
246                        !
247                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn
248                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta
249                     END DO
250                  END DO
251               END DO
252            END DO
253         ENDIF
254         !
255      ENDIF
256      !
257   END SUBROUTINE tra_nxt_fix
258
259
260   SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt )
261      !!----------------------------------------------------------------------
262      !!                   ***  ROUTINE tra_nxt_vvl  ***
263      !!
264      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
265      !!                and swap the tracer fields.
266      !!
267      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
268      !!              - save in (ta,sa) a thickness weighted average over the three
269      !!             time levels which will be used to compute rdn and thus the semi-
270      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
271      !!              - swap tracer fields to prepare the next time_step.
272      !!                This can be summurized for tempearture as:
273      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
274      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
275      !!             ztm = 0                                                       otherwise
276      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
277      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
278      !!             tn  = ta
279      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
280      !!
281      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
282      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
283      !!----------------------------------------------------------------------
284      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
285      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
286      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
287      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
288      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
289      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
290      !!     
291      INTEGER  ::   ji, jj, jk, jn                 ! dummy loop indices
292      REAL(wp) ::   ztc_a , ztc_n , ztc_b          ! temporary scalar
293      REAL(wp) ::   ztc_f , ztc_d , ztc_m          !    -         -
294      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a         !    -         -
295      REAL(wp) ::   ze3t_f, ze3t_d, ze3t_m         !    -         -
296      REAL     ::   zfact1, zfact2                 !    -         -
297      !!----------------------------------------------------------------------
298
299      IF( kt == nit000 ) THEN
300         IF(lwp) WRITE(numout,*)
301         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
302         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
303      ENDIF
304      !
305      !
306      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step
307         !                                           ! (only swap)
308         DO jn = 1, kjpt
309            DO jk = 1, jpkm1
310               ptn(:,:,jk,jn) = pta(:,:,jk,jn)       ! ptb <-- ptn
311            END DO
312         END DO
313         !                                             
314      ELSE                                           ! general case (Leapfrog + Asselin filter)
315         !
316         !                                           ! ----------------------- !
317         IF( ln_dynhpg_imp ) THEN                    ! semi-implicite hpg case !
318            !                                        ! ----------------------- !
319            DO jn = 1, kjpt                         
320               DO jk = 1, jpkm1
321                  DO jj = 1, jpj
322                     DO ji = 1, jpi
323                        ze3t_b = fse3t_b(ji,jj,jk)
324                        ze3t_n = fse3t_n(ji,jj,jk)
325                        ze3t_a = fse3t_a(ji,jj,jk)
326                        !                                         ! tracer content at Before, now and after
327                        ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
328                        ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n 
329                        ztc_a  = pta(ji,jj,jk,jn) * ze3t_a 
330                        !                                         ! Time laplacian on tracer contents
331                        !                                         ! used for both Asselin and Brown & Campana filters
332                        ze3t_d =  ze3t_b - 2.* ze3t_n + ze3t_a 
333                        ztc_d  =  ztc_b  - 2.* ztc_n  + ztc_a 
334                        !                                         ! Asselin Filter on thicknesses and tracer contents
335                        ztc_f  = ztc_n  + atfp * ztc_d
336                        ztc_m  = ztc_n  + rbcp * ztc_d
337                        !                                         
338                        ze3t_f = 1.0 / ( ze3t_n + atfp * ze3t_d )
339                        ze3t_m = 1.0 / ( ze3t_n + rbcp * ze3t_d )
340                        !                                         ! swap of arrays
341                        ptb(ji,jj,jk,jn) = ztc_f * ze3t_f              ! ptb <-- ptn filtered
342                        pta(ji,jj,jk,jn) = ztc_m * ze3t_m              ! pta <-- Brown & Campana average
343                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)            ! ptn <-- pta
344                     END DO
345                  END DO
346               END DO
347            END DO
348            !                                           ! ----------------------- !
349         ELSE                                           !    explicit hpg case    !
350            !                                           ! ----------------------- !
351            IF( cdtype == 'TRA' ) THEN
352               !             
353               DO jn = 1, kjpt     
354                  DO jk = 1, jpkm1
355                     zfact1 = atfp * rdttra(jk)
356                     zfact2 = zfact1 / rau0
357                     DO jj = 1, jpj
358                        DO ji = 1, jpi
359                           ze3t_b = fse3t_b(ji,jj,jk)
360                           ze3t_n = fse3t_n(ji,jj,jk)
361                           ze3t_a = fse3t_a(ji,jj,jk)
362                           !                                         ! tracer content at Before, now and after
363                           ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
364                           ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
365                           ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
366                           !
367                           ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a )
368                           ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b )
369
370                           IF( jk == 1 ) THEN
371                               ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
372                               IF( jn == jp_tem )   ztc_f  = ztc_f  - zfact1 * ( sbc_hc_n(ji,jj) - sbc_hc_b(ji,jj) )
373                               IF( jn == jp_sal )   ztc_f  = ztc_f  - zfact1 * ( sbc_sc_n(ji,jj) - sbc_sc_b(ji,jj) )
374                           ENDIF
375                           IF( jn == jp_tem .AND. ln_traqsr .AND. jk <= nksr )  &
376                           &                        ztc_f  = ztc_f  - zfact1 * ( qsr_hc_n(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
377
378                           ze3t_f = 1.e0 / ze3t_f
379                           ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered
380                           ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta
381                        END DO
382                     END DO
383                  END DO
384               END DO
385               !
386            ELSE IF( cdtype == 'TRC' ) THEN
387               !             
388               DO jn = 1, kjpt     
389                  DO jk = 1, jpkm1
390                     DO jj = 1, jpj
391                        DO ji = 1, jpi
392                           ze3t_b = fse3t_b(ji,jj,jk)
393                           ze3t_n = fse3t_n(ji,jj,jk)
394                           ze3t_a = fse3t_a(ji,jj,jk)
395                           !                                         ! tracer content at Before, now and after
396                           ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
397                           ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
398                           ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
399                           !
400                           ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a )
401                           ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b  )
402
403                           ze3t_f = 1.e0 / ze3t_f
404                           ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered
405                           ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta
406                        END DO
407                     END DO
408                  END DO
409               END DO
410               !
411            END IF
412            !
413         ENDIF
414         !
415      ENDIF
416      !
417   END SUBROUTINE tra_nxt_vvl
418
419   !!======================================================================
420END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.