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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 18.1 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 sbcrnf          ! river runoffs
30   USE sbcisf          ! ice shelf melting
31   USE zdf_oce         ! ocean vertical mixing
32   USE domvvl          ! variable volume
33   USE trd_oce         ! trends: ocean variables
34   USE trdtra          ! trends manager: tracers
35   USE traqsr          ! penetrative solar radiation (needed for nksr)
36   USE phycst          ! physical constant
37   USE ldftra          ! lateral physics on tracers
38   USE ldfslp
39   USE bdy_oce   , ONLY: ln_bdy
40   USE bdytra          ! open boundary condition (bdy_tra routine)
41   !
42   USE in_out_manager  ! I/O manager
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE prtctl          ! Print control
45   USE wrk_nemo        ! Memory allocation
46   USE timing          ! Timing
47#if defined key_agrif
48   USE agrif_opa_interp
49#endif
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC   tra_nxt       ! routine called by step.F90
55   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
56   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
57
58   !! * Substitutions
59#  include "vectopt_loop_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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 one-way open boundaries (ln_bdy=T),
82      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
83      !!
84      !!              - Update lateral boundary conditions on AGRIF children
85      !!             domains (lk_agrif=T)
86      !!
87      !! ** Action  : - tsb & tsn ready for the next time step
88      !!----------------------------------------------------------------------
89      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
90      !!
91      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
92      REAL(wp) ::   zfact            ! local scalars
93      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
94      !!----------------------------------------------------------------------
95      !
96      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt')
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      ENDIF
103
104      ! Update after tracer on domain lateral boundaries
105      !
106#if defined key_agrif
107      CALL Agrif_tra                     ! AGRIF zoom boundaries
108#endif
109      !
110      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign)
111      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
112      !
113      IF( ln_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
114 
115      ! set time step size (Euler/Leapfrog)
116      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =     rdt      ! at nit000             (Euler)
117      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
118      ENDIF
119
120      ! trends computation initialisation
121      IF( l_trdtra )   THEN                   
122         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
123!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
124         DO jk = 1, jpk
125            DO jj = 1, jpj
126               DO ji = 1, jpi
127                  ztrdt(ji,jj,jk) = 0._wp 
128                  ztrds(ji,jj,jk) = 0._wp
129               END DO
130            END DO
131         END DO
132         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
133            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
134            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
135         ENDIF
136         ! total trend for the non-time-filtered variables.
137            zfact = 1.0 / rdt
138!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
139         DO jk = 1, jpkm1
140            DO jj = 1, jpj
141               DO ji = 1, jpi
142                  ztrdt(ji,jj,jk) = ( tsa(ji,jj,jk,jp_tem) - tsn(ji,jj,jk,jp_tem) ) * zfact
143                  ztrds(ji,jj,jk) = ( tsa(ji,jj,jk,jp_sal) - tsn(ji,jj,jk,jp_sal) ) * zfact
144               END DO
145            END DO
146         END DO
147         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
148         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
149         ! Store now fields before applying the Asselin filter
150         ! in order to calculate Asselin filter trend later.
151!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
152         DO jk = 1, jpkm1
153            DO jj = 1, jpj
154               DO ji = 1, jpi
155                  ztrdt(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
156                  ztrds(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
157               END DO
158            END DO
159         END DO
160      ENDIF
161
162      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
163         DO jn = 1, jpts
164!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
165            DO jk = 1, jpkm1
166               DO jj = 1, jpj
167                  DO ji = 1, jpi
168                     tsn(ji,jj,jk,jn) = tsa(ji,jj,jk,jn)   
169                  END DO
170               END DO
171            END DO
172         END DO
173         !
174      ELSE                                            ! Leap-Frog + Asselin filter time stepping
175         !
176         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface
177         ELSE                   ;   CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa,   &
178           &                                                                sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
179         ENDIF
180         !
181         DO jn = 1, jpts
182            CALL lbc_lnk( tsb(:,:,:,jn), 'T', 1._wp ) 
183            CALL lbc_lnk( tsn(:,:,:,jn), 'T', 1._wp )
184            CALL lbc_lnk( tsa(:,:,:,jn), 'T', 1._wp )
185         END DO
186      ENDIF     
187      !
188      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
189!$OMP PARALLEL DO schedule(static) private(jk, zfact)
190         DO jk = 1, jpkm1
191            DO jj = 1, jpj
192               DO ji = 1, jpi
193                  zfact = 1._wp / r2dt             
194                  ztrdt(ji,jj,jk) = ( tsb(ji,jj,jk,jp_tem) - ztrdt(ji,jj,jk) ) * zfact
195                  ztrds(ji,jj,jk) = ( tsb(ji,jj,jk,jp_sal) - ztrds(ji,jj,jk) ) * zfact
196               END DO
197            END DO
198         END DO
199         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
200         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
201         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
202      END IF
203      !
204      !                        ! control print
205      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
206         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
207      !
208      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
209      !
210   END SUBROUTINE tra_nxt
211
212
213   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
214      !!----------------------------------------------------------------------
215      !!                   ***  ROUTINE tra_nxt_fix  ***
216      !!
217      !! ** Purpose :   fixed volume: apply the Asselin time filter and
218      !!                swap the tracer fields.
219      !!
220      !! ** Method  : - Apply a Asselin time filter on now fields.
221      !!              - swap tracer fields to prepare the next time_step.
222      !!
223      !! ** Action  : - tsb & tsn ready for the next time step
224      !!----------------------------------------------------------------------
225      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
226      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
227      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
228      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
229      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
230      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
231      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
232      !
233      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
234      REAL(wp) ::   ztn, ztd         ! local scalars
235      !!----------------------------------------------------------------------
236      !
237      IF( kt == kit000 )  THEN
238         IF(lwp) WRITE(numout,*)
239         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
240         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
241      ENDIF
242      !
243      DO jn = 1, kjpt
244         !
245!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,ztn,ztd)
246         DO jk = 1, jpkm1
247            DO jj = 2, jpjm1
248               DO ji = fs_2, fs_jpim1
249                  ztn = ptn(ji,jj,jk,jn)                                   
250                  ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn)  ! time laplacian on tracers
251                  !
252                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                      ! ptb <-- filtered ptn
253                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                      ! ptn <-- pta
254               END DO
255           END DO
256         END DO
257         !
258      END DO
259      !
260   END SUBROUTINE tra_nxt_fix
261
262
263   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, 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      !!              - swap tracer fields to prepare the next time_step.
272      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
273      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
274      !!             tn  = ta
275      !!
276      !! ** Action  : - tsb & tsn ready for the next time step
277      !!----------------------------------------------------------------------
278      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
279      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
280      REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step
281      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
282      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
283      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
284      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
285      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
286      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc   ! surface tracer content
287      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
288      !
289      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
290      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
291      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
292      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
293      !!----------------------------------------------------------------------
294      !
295      IF( kt == kit000 )  THEN
296         IF(lwp) WRITE(numout,*)
297         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
298         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
299      ENDIF
300      !
301      IF( cdtype == 'TRA' )  THEN   
302         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
303         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
304         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
305      ELSE                          ! passive tracers case
306         ll_traqsr  = .FALSE.          ! NO solar penetration
307         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ? 
308         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??
309      ENDIF
310      !
311      DO jn = 1, kjpt     
312!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zfact1,zfact2,ze3t_b,ze3t_n,ze3t_a,ze3t_d,ze3t_f,ztc_b,ztc_n,ztc_a,ztc_d,ztc_f)
313         DO jk = 1, jpkm1
314            zfact1 = atfp * p2dt
315            zfact2 = zfact1 * r1_rau0
316            DO jj = 2, jpjm1
317               DO ji = fs_2, fs_jpim1
318                  ze3t_b = e3t_b(ji,jj,jk)
319                  ze3t_n = e3t_n(ji,jj,jk)
320                  ze3t_a = e3t_a(ji,jj,jk)
321                  !                                         ! tracer content at Before, now and after
322                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
323                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
324                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
325                  !
326                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
327                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
328                  !
329                  ze3t_f = ze3t_n + atfp * ze3t_d
330                  ztc_f  = ztc_n  + atfp * ztc_d
331                  !
332                  IF( jk == mikt(ji,jj) ) THEN           ! first level
333                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
334                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
335                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
336                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
337                  ENDIF
338                  !
339                  ! solar penetration (temperature only)
340                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
341                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
342                     !
343                  ! river runoff
344                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
345                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
346                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
347                     !
348                  ! ice shelf
349                  IF( ll_isf ) THEN
350                     ! level fully include in the Losch_2008 ice shelf boundary layer
351                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
352                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
353                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
354                     ! level partially include in Losch_2008 ice shelf boundary layer
355                     IF ( jk == misfkb(ji,jj) )                                                   &
356                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
357                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
358                  END IF
359                  !
360                  ze3t_f = 1.e0 / ze3t_f
361                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
362                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
363                  !
364               END DO
365            END DO
366         END DO
367         !
368      END DO
369      !
370   END SUBROUTINE tra_nxt_vvl
371
372   !!======================================================================
373END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.