source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

  • Property svn:keywords set to Id
File size: 18.0 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 zdf_oce         ! ocean vertical mixing
31   USE domvvl          ! variable volume
32   USE dynspg_oce      ! surface     pressure gradient variables
33   USE dynhpg          ! hydrostatic pressure gradient
34   USE trd_oce         ! trends: ocean variables
35   USE trdtra          ! trends manager: tracers
36   USE traqsr          ! penetrative solar radiation (needed for nksr)
37   USE phycst          ! physical constant
38   USE ldftra_oce      ! lateral physics on tracers
39   USE bdy_oce         ! BDY open boundary condition variables
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_update
49   USE agrif_opa_interp
50#endif
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_bdy=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_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - 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._wp )      ! local domain boundaries  (T-point, unchanged sign)
113      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
114      !
115#if defined key_bdy 
116      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
117#endif
118#if defined key_agrif
119      CALL Agrif_tra                     ! AGRIF zoom boundaries
120#endif
121 
122      ! set time step size (Euler/Leapfrog)
123      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
124      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
125      ENDIF
126
127      ! trends computation initialisation
128      IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter
129         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
130         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
131         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
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      ENDIF
137
138      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
139         DO jn = 1, jpts
140            DO jk = 1, jpkm1
141               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
142            END DO
143         END DO
144      ELSE                                            ! Leap-Frog + Asselin filter time stepping
145         !
146         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   &
147           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)
148         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
149         ENDIF
150      ENDIF 
151      !
152#if defined key_agrif
153      ! Update tracer at AGRIF zoom boundaries
154      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
155#endif     
156      !
157      ! trends computation
158      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
159         DO jk = 1, jpkm1
160            zfact = 1._wp / r2dtra(jk)             
161            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
162            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
163         END DO
164         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
165         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
166         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
167      END IF
168      !
169      !                        ! control print
170      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
171         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
172      !
173      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
174      !
175   END SUBROUTINE tra_nxt
176
177
178   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
179      !!----------------------------------------------------------------------
180      !!                   ***  ROUTINE tra_nxt_fix  ***
181      !!
182      !! ** Purpose :   fixed volume: apply the Asselin time filter and
183      !!                swap the tracer fields.
184      !!
185      !! ** Method  : - Apply a Asselin time filter on now fields.
186      !!              - save in (ta,sa) an average over the three time levels
187      !!             which will be used to compute rdn and thus the semi-implicit
188      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
189      !!              - swap tracer fields to prepare the next time_step.
190      !!                This can be summurized for tempearture as:
191      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
192      !!             ztm = 0                                   otherwise
193      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
194      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
195      !!             tn  = ta 
196      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
197      !!
198      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
199      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
200      !!----------------------------------------------------------------------
201      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
202      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
203      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
204      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
205      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
206      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
207      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
208      !
209      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
210      LOGICAL  ::   ll_tra_hpg       ! local logical
211      REAL(wp) ::   ztn, ztd         ! local scalars
212      !!----------------------------------------------------------------------
213
214      IF( kt == kit000 )  THEN
215         IF(lwp) WRITE(numout,*)
216         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
217         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
218      ENDIF
219      !
220      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
221      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
222      ENDIF
223      !
224      DO jn = 1, kjpt
225         !
226         DO jk = 1, jpkm1
227            DO jj = 1, jpj
228               DO ji = 1, jpi
229                  ztn = ptn(ji,jj,jk,jn)                                   
230                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
231                  !
232                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
233                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
234                  !
235                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
236               END DO
237           END DO
238         END DO
239         !
240      END DO
241      !
242   END SUBROUTINE tra_nxt_fix
243
244
245   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
246      !!----------------------------------------------------------------------
247      !!                   ***  ROUTINE tra_nxt_vvl  ***
248      !!
249      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
250      !!                and swap the tracer fields.
251      !!
252      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
253      !!              - save in (ta,sa) a thickness weighted average over the three
254      !!             time levels which will be used to compute rdn and thus the semi-
255      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
256      !!              - swap tracer fields to prepare the next time_step.
257      !!                This can be summurized for tempearture as:
258      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
259      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
260      !!             ztm = 0                                                       otherwise
261      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
262      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
263      !!             tn  = ta
264      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
265      !!
266      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
267      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
268      !!----------------------------------------------------------------------
269      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
270      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
271      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
272      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
273      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
274      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
275      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
276      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
277      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
278      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
279
280      !!     
281      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical
282      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
283      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
284      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
285      !!----------------------------------------------------------------------
286      !
287      IF( kt == kit000 )  THEN
288         IF(lwp) WRITE(numout,*)
289         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
290         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
291      ENDIF
292      !
293      IF( cdtype == 'TRA' )  THEN   
294         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
295         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
296         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
297      ELSE                         
298         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
299         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
300         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
301      ENDIF
302      !
303      DO jn = 1, kjpt     
304         DO jk = 1, jpkm1
305            zfact1 = atfp * p2dt(jk)
306            zfact2 = zfact1 / rau0
307            DO jj = 1, jpj
308               DO ji = 1, jpi
309                  ze3t_b = fse3t_b(ji,jj,jk)
310                  ze3t_n = fse3t_n(ji,jj,jk)
311                  ze3t_a = fse3t_a(ji,jj,jk)
312                  !                                         ! tracer content at Before, now and after
313                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
314                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
315                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
316                  !
317                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
318                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
319                  !
320                  ze3t_f = ze3t_n + atfp * ze3t_d
321                  ztc_f  = ztc_n  + atfp * ztc_d
322                  !
323                  IF( jk == 1 ) THEN           ! first level
324                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) )
325                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
326                  ENDIF
327
328                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
329                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
330
331                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )   &            ! river runoffs
332                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
333                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
334
335                  ze3t_f = 1.e0 / ze3t_f
336                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
337                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
338                  !
339                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
340                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
341                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
342                  ENDIF
343               END DO
344            END DO
345         END DO
346         !
347      END DO
348      !
349   END SUBROUTINE tra_nxt_vvl
350
351   !!======================================================================
352END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.