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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago

Ediag branche: #927 split TRA/DYN trd computation
  • Property svn:keywords set to Id
File size: 17.4 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 phycst          ! physical constant
32   USE ldftra_oce      ! ocean active tracers lateral physics
33   USE dynspg_oce      ! surface     pressure gradient variables
34   USE dynhpg          ! hydrostatic pressure gradient
35   USE traqsr          ! penetrative solar radiation (needed for nksr)
36   USE trd_oce         ! trends: ocean variables
37   USE trdtra          ! trends manager: tracers
38   USE obc_oce         ! open boundary condition variable
39   USE obctra          ! open boundary condition (obc_tra routine)
40   USE bdy_oce
41   USE bdytra          ! open boundary condition (bdy_tra routine)
42   USE in_out_manager  ! I/O manager
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE prtctl          ! Print control
45#if defined key_agrif
46   USE agrif_opa_update
47   USE agrif_opa_interp
48#endif
49   USE wrk_nemo        ! Memory allocation
50   USE timing          ! 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   !! * Substitutions
62#  include "domzgr_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
65   !! $Id$
66   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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 one-way open boundaries (lk_obc=T),
85      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
86      !!
87      !!              - Update lateral boundary conditions on AGRIF children
88      !!             domains (lk_agrif=T)
89      !!
90      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
91      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
92      !!----------------------------------------------------------------------
93      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
94      !!
95      INTEGER  ::   jk, jn    ! dummy loop indices
96      REAL(wp) ::   zfact     ! local scalars
97      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
98      !!----------------------------------------------------------------------
99      !
100      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt')
101      !
102      IF( kt == nit000 ) THEN
103         IF(lwp) WRITE(numout,*)
104         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
105         IF(lwp) WRITE(numout,*) '~~~~~~~'
106         !
107         rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)      ! Brown & Campana parameter for semi-implicit hpg
108      ENDIF
109
110      ! Update after tracer on domain lateral boundaries
111      !
112      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )      ! local domain boundaries  (T-point, unchanged sign)
113      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
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      ! set time step size (Euler/Leapfrog)
126      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
127      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
128      ENDIF
129
130      ! trends computation initialisation
131      IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter
132         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
133         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
134         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
135         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
136            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
137            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
138         ENDIF
139      ENDIF
140
141      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
142         DO jn = 1, jpts
143            DO jk = 1, jpkm1
144               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
145            END DO
146         END DO
147      ELSE                                            ! Leap-Frog + Asselin filter time stepping
148         !
149         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
150         ELSE                 ;   CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
151         ENDIF
152      ENDIF 
153      !
154#if defined key_agrif
155      ! Update tracer at AGRIF zoom boundaries
156      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
157#endif     
158      !
159      ! trends computation
160      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
161         DO jk = 1, jpkm1
162            zfact = 1.e0 / r2dtra(jk)             
163            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
164            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
165         END DO
166         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
167         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
168         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
169      END IF
170      !
171      !                        ! control print
172      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
173         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
174      !
175      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
176      !
177   END SUBROUTINE tra_nxt
178
179
180   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
181      !!----------------------------------------------------------------------
182      !!                   ***  ROUTINE tra_nxt_fix  ***
183      !!
184      !! ** Purpose :   fixed volume: apply the Asselin time filter and
185      !!                swap the tracer fields.
186      !!
187      !! ** Method  : - Apply a Asselin time filter on now fields.
188      !!              - save in (ta,sa) an average over the three time levels
189      !!             which will be used to compute rdn and thus the semi-implicit
190      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
191      !!              - swap tracer fields to prepare the next time_step.
192      !!                This can be summurized for tempearture as:
193      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
194      !!             ztm = 0                                   otherwise
195      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
196      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
197      !!             tn  = ta 
198      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
199      !!
200      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
201      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
202      !!----------------------------------------------------------------------
203      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
204      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
205      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
206      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
207      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
208      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
209      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
210      !
211      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
212      LOGICAL  ::   ll_tra_hpg       ! local logical
213      REAL(wp) ::   ztn, ztd         ! local scalars
214      !!----------------------------------------------------------------------
215
216      IF( kt == kit000 )  THEN
217         IF(lwp) WRITE(numout,*)
218         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
219         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
220      ENDIF
221      !
222      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
223      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
224      ENDIF
225      !
226      DO jn = 1, kjpt
227         !
228         DO jk = 1, jpkm1
229            DO jj = 1, jpj
230               DO ji = 1, jpi
231                  ztn = ptn(ji,jj,jk,jn)                                   
232                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
233                  !
234                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
235                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
236                  !
237                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
238               END DO
239           END DO
240         END DO
241         !
242      END DO
243      !
244   END SUBROUTINE tra_nxt_fix
245
246
247   SUBROUTINE tra_nxt_vvl( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
248      !!----------------------------------------------------------------------
249      !!                   ***  ROUTINE tra_nxt_vvl  ***
250      !!
251      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
252      !!                and swap the tracer fields.
253      !!
254      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
255      !!              - save in (ta,sa) a thickness weighted average over the three
256      !!             time levels which will be used to compute rdn and thus the semi-
257      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
258      !!              - swap tracer fields to prepare the next time_step.
259      !!                This can be summurized for tempearture as:
260      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
261      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
262      !!             ztm = 0                                                       otherwise
263      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
264      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
265      !!             tn  = ta
266      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
267      !!
268      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
269      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
270      !!----------------------------------------------------------------------
271      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
272      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
273      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
274      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
275      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
276      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
277      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
278      !!     
279      LOGICAL  ::   ll_tra, ll_tra_hpg, ll_traqsr   ! local logical
280      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
281      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
282      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
283      !!----------------------------------------------------------------------
284      !
285      IF( kt == kit000 )  THEN
286         IF(lwp) WRITE(numout,*)
287         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
288         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
289      ENDIF
290      !
291      IF( cdtype == 'TRA' )  THEN   
292         ll_tra     = .TRUE.           ! active tracers case 
293         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
294         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
295      ELSE                         
296         ll_tra     = .FALSE.          ! passive tracers case
297         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
298         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
299      ENDIF
300      !
301      DO jn = 1, kjpt     
302         DO jk = 1, jpkm1
303            zfact1 = atfp * rdttra(jk)
304            zfact2 = zfact1 / rau0
305            DO jj = 1, jpj
306               DO ji = 1, jpi
307                  ze3t_b = fse3t_b(ji,jj,jk)
308                  ze3t_n = fse3t_n(ji,jj,jk)
309                  ze3t_a = fse3t_a(ji,jj,jk)
310                  !                                         ! tracer content at Before, now and after
311                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
312                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
313                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
314                  !
315                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
316                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
317                  !
318                  ze3t_f = ze3t_n + atfp * ze3t_d
319                  ztc_f  = ztc_n  + atfp * ztc_d
320                  !
321                  IF( ll_tra .AND. jk == 1 ) THEN           ! first level only for T & S
322                      ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
323                      ztc_f  = ztc_f  - zfact1 * ( sbc_tsc(ji,jj,jn) - sbc_tsc_b(ji,jj,jn) )
324                  ENDIF
325                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
326                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
327
328                   ze3t_f = 1.e0 / ze3t_f
329                   ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
330                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
331                   !
332                   IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
333                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
334                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
335                   ENDIF
336               END DO
337            END DO
338         END DO
339         !
340      END DO
341      !
342   END SUBROUTINE tra_nxt_vvl
343
344   !!======================================================================
345END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.