source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 11134

Last change on this file since 11134 was 11134, checked in by jcastill, 18 months ago

Full set of changes as in the original branch

File size: 20.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 sbcrnf          ! river runoffs
30   USE sbcisf          ! ice shelf melting/freezing
31   USE zdf_oce         ! ocean vertical mixing
32   USE domvvl          ! variable volume
33   USE dynspg_oce      ! surface     pressure gradient variables
34   USE dynhpg          ! hydrostatic pressure gradient
35   USE trd_oce         ! trends: ocean variables
36   USE trdtra          ! trends manager: tracers
37   USE traqsr          ! penetrative solar radiation (needed for nksr)
38   USE phycst          ! physical constant
39   USE ldftra_oce      ! lateral physics on tracers
40   USE bdy_oce         ! BDY open boundary condition variables
41   USE bdytra          ! open boundary condition (bdy_tra routine)
42   !
43   USE in_out_manager  ! I/O manager
44   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
45   USE prtctl          ! Print control
46   USE wrk_nemo        ! Memory allocation
47   USE timing          ! Timing
48#if defined key_agrif
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   INTEGER  ::   warn_1, warn_2   ! indicators for warning statement
61
62   !! * Substitutions
63#  include "domzgr_substitute.h90"
64   !!----------------------------------------------------------------------
65   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
66   !! $Id$
67   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
68   !!----------------------------------------------------------------------
69CONTAINS
70
71   SUBROUTINE tra_nxt( kt )
72      !!----------------------------------------------------------------------
73      !!                   ***  ROUTINE tranxt  ***
74      !!
75      !! ** Purpose :   Apply the boundary condition on the after temperature 
76      !!             and salinity fields, achieved the time stepping by adding
77      !!             the Asselin filter on now fields and swapping the fields.
78      !!
79      !! ** Method  :   At this stage of the computation, ta and sa are the
80      !!             after temperature and salinity as the time stepping has
81      !!             been performed in trazdf_imp or trazdf_exp module.
82      !!
83      !!              - Apply lateral boundary conditions on (ta,sa)
84      !!             at the local domain   boundaries through lbc_lnk call,
85      !!             at the one-way open boundaries (lk_bdy=T),
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, jn, ji, jj     ! dummy loop indices
97      REAL(wp) ::   zfact, zfreeze     ! local scalars
98      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
99      !!----------------------------------------------------------------------
100      !
101      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt')
102      !
103      IF( kt == nit000 ) THEN
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
106         IF(lwp) WRITE(numout,*) '~~~~~~~'
107         !
108         rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp)      ! Brown & Campana parameter for semi-implicit hpg
109      ENDIF
110
111      ! Update after tracer on domain lateral boundaries
112      !
113#if defined key_agrif
114      CALL Agrif_tra                     ! AGRIF zoom boundaries
115#endif
116      !
117      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign)
118      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
119      !
120#if defined key_bdy 
121      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
122#endif
123 
124      ! set time step size (Euler/Leapfrog)
125      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
126      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
127      ENDIF
128
129#if ( ! defined key_lim3 && ! defined key_lim2 && ! key_cice )
130      IF ( kt == nit000 ) warn_1=0
131      warn_2=0
132      DO jk = 1, jpkm1
133         DO jj = 1, jpj
134            DO ji = 1, jpi
135               IF ( tsa(ji,jj,jk,jp_tem) .lt. 0.0 ) THEN
136                  ! calculate freezing point
137                  zfreeze = ( -0.0575_wp + 1.710523E-3 * Sqrt(Abs(tsn(ji,jj,jk,jp_sal)))   & 
138                            - 2.154996E-4 * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) - 7.53E-4 * ( 10.0_wp + fsdept(ji,jj,jk) )
139                  IF ( tsa(ji,jj,jk,jp_tem) .lt. zfreeze ) THEN
140                     tsa(ji,jj,jk,jp_tem)=zfreeze
141                     warn_2=1
142                  ENDIF
143               ENDIF
144            END DO
145         END DO
146      END DO
147      CALL mpp_max(warn_1)
148      CALL mpp_max(warn_2)
149      IF ( (warn_1 == 0) .and. (warn_2 /= 0) ) THEN
150         IF(lwp) THEN
151            CALL ctl_warn( ' Temperatures dropping below freezing point, ', &
152                      &    ' being forced to freezing point, no longer conservative' ) 
153         ENDIF
154         warn_1=1
155      ENDIF
156#endif
157
158      ! trends computation initialisation
159      IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter
160         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
161         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
162         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
163         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
164            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
165            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
166         ENDIF
167      ENDIF
168
169      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
170         DO jn = 1, jpts
171            DO jk = 1, jpkm1
172               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
173            END DO
174         END DO
175      ELSE                                            ! Leap-Frog + Asselin filter time stepping
176         !
177         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   &
178           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)
179         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
180         ENDIF
181      ENDIF     
182      !
183     ! trends computation
184      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
185         DO jk = 1, jpkm1
186            zfact = 1._wp / r2dtra(jk)             
187            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
188            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
189         END DO
190         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
191         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
192         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
193      END IF
194      !
195      !                        ! control print
196      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
197         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
198      !
199      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
200      !
201   END SUBROUTINE tra_nxt
202
203
204   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
205      !!----------------------------------------------------------------------
206      !!                   ***  ROUTINE tra_nxt_fix  ***
207      !!
208      !! ** Purpose :   fixed volume: apply the Asselin time filter and
209      !!                swap the tracer fields.
210      !!
211      !! ** Method  : - Apply a Asselin time filter on now fields.
212      !!              - save in (ta,sa) an average over the three time levels
213      !!             which will be used to compute rdn and thus the semi-implicit
214      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
215      !!              - swap tracer fields to prepare the next time_step.
216      !!                This can be summurized for tempearture as:
217      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
218      !!             ztm = 0                                   otherwise
219      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
220      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
221      !!             tn  = ta 
222      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
223      !!
224      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
225      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
226      !!----------------------------------------------------------------------
227      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
228      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
229      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
230      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
231      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
232      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
233      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
234      !
235      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
236      LOGICAL  ::   ll_tra_hpg       ! local logical
237      REAL(wp) ::   ztn, ztd         ! local scalars
238      !!----------------------------------------------------------------------
239
240      IF( kt == kit000 )  THEN
241         IF(lwp) WRITE(numout,*)
242         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
243         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
244      ENDIF
245      !
246      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg   
247      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
248      ENDIF
249      !
250      DO jn = 1, kjpt
251         !
252         DO jk = 1, jpkm1
253            DO jj = 1, jpj
254               DO ji = 1, jpi
255                  ztn = ptn(ji,jj,jk,jn)                                   
256                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
257                  !
258                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn
259                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
260                  !
261                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
262               END DO
263           END DO
264         END DO
265         !
266      END DO
267      !
268   END SUBROUTINE tra_nxt_fix
269
270
271   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
272      !!----------------------------------------------------------------------
273      !!                   ***  ROUTINE tra_nxt_vvl  ***
274      !!
275      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
276      !!                and swap the tracer fields.
277      !!
278      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
279      !!              - save in (ta,sa) a thickness weighted average over the three
280      !!             time levels which will be used to compute rdn and thus the semi-
281      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
282      !!              - swap tracer fields to prepare the next time_step.
283      !!                This can be summurized for tempearture as:
284      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
285      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
286      !!             ztm = 0                                                       otherwise
287      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
288      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
289      !!             tn  = ta
290      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
291      !!
292      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
293      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
294      !!----------------------------------------------------------------------
295      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
296      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
297      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
298      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
299      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
300      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
301      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
302      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
303      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
304      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
305
306      !!     
307      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf, ll_isf   ! local logical
308      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
309      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
310      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
311      !!----------------------------------------------------------------------
312      !
313      IF( kt == kit000 )  THEN
314         IF(lwp) WRITE(numout,*)
315         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
316         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
317      ENDIF
318      !
319      IF( cdtype == 'TRA' )  THEN   
320         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
321         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
322         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
323         IF (nn_isf .GE. 1) THEN
324            ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing
325         ELSE
326            ll_isf = .FALSE.
327         END IF
328      ELSE                         
329         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
330         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
331         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
332         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing
333      ENDIF
334      !
335      DO jn = 1, kjpt     
336         DO jk = 1, jpkm1
337            zfact1 = atfp * p2dt(jk)
338            zfact2 = zfact1 / rau0
339            DO jj = 1, jpj
340               DO ji = 1, jpi
341                  ze3t_b = fse3t_b(ji,jj,jk)
342                  ze3t_n = fse3t_n(ji,jj,jk)
343                  ze3t_a = fse3t_a(ji,jj,jk)
344                  !                                         ! tracer content at Before, now and after
345                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
346                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
347                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
348                  !
349                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
350                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
351                  !
352                  ze3t_f = ze3t_n + atfp * ze3t_d
353                  ztc_f  = ztc_n  + atfp * ztc_d
354                  !
355                  IF( jk == mikt(ji,jj) ) THEN           ! first level
356                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
357                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
358                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
359                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
360                  ENDIF
361
362                  ! solar penetration (temperature only)
363                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
364                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
365
366                  ! river runoff
367                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
368                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
369                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
370
371                  ! ice shelf
372                  IF( ll_isf ) THEN
373                     ! level fully include in the Losch_2008 ice shelf boundary layer
374                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
375                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
376                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
377                     ! level partially include in Losch_2008 ice shelf boundary layer
378                     IF ( jk == misfkb(ji,jj) )                                                   &
379                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
380                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
381                  END IF
382
383                  ze3t_f = 1.e0 / ze3t_f
384                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
385                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
386                  !
387                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
388                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
389                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
390                  ENDIF
391               END DO
392            END DO
393         END DO
394         !
395      END DO
396      !
397   END SUBROUTINE tra_nxt_vvl
398
399   !!======================================================================
400END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.