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

source: branches/UKMO/dev_r5518_amm15_test/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

Last change on this file was 6662, checked in by jgraham, 8 years ago

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