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

Last change on this file since 2236 was 2236, checked in by cetlod, 13 years ago

First guess of NEMO_v3.3

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