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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 4451

Last change on this file since 4451 was 4451, checked in by trackstand2, 10 years ago

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