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

source: branches/2011/dev_UKM0_2011/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 3062

Last change on this file since 3062 was 3062, checked in by rfurner, 12 years ago

ticket #885. added in changes from branches/2011/UKMO_MERCATOR_obc_bdy_merge@2888
  • Property svn:keywords set to Id
File size: 16.8 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_oce
39   USE bdytra          ! open boundary condition (bdy_tra 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#if defined key_agrif
46   USE agrif_opa_update
47   USE agrif_opa_interp
48#endif
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   tra_nxt       ! routine called by step.F90
54   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
55   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
56
57   REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg
58
59   !! * Substitutions
60#  include "domzgr_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE tra_nxt( kt )
69      !!----------------------------------------------------------------------
70      !!                   ***  ROUTINE tranxt  ***
71      !!
72      !! ** Purpose :   Apply the boundary condition on the after temperature 
73      !!             and salinity fields, achieved the time stepping by adding
74      !!             the Asselin filter on now fields and swapping the fields.
75      !!
76      !! ** Method  :   At this stage of the computation, ta and sa are the
77      !!             after temperature and salinity as the time stepping has
78      !!             been performed in trazdf_imp or trazdf_exp module.
79      !!
80      !!              - Apply lateral boundary conditions on (ta,sa)
81      !!             at the local domain   boundaries through lbc_lnk call,
82      !!             at the one-way open boundaries (lk_obc=T),
83      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
84      !!
85      !!              - Update lateral boundary conditions on AGRIF children
86      !!             domains (lk_agrif=T)
87      !!
88      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
89      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
90      !!----------------------------------------------------------------------
91      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
92      !!
93      INTEGER  ::   jk, jn    ! dummy loop indices
94      REAL(wp) ::   zfact     ! local scalars
95      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
96      !!----------------------------------------------------------------------
97
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
101         IF(lwp) WRITE(numout,*) '~~~~~~~'
102         !
103         rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)      ! Brown & Campana parameter for semi-implicit hpg
104      ENDIF
105
106      ! Update after tracer on domain lateral boundaries
107      !
108      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )      ! local domain boundaries  (T-point, unchanged sign)
109      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
110      !
111#if defined key_obc || defined key_bdy || defined key_agrif
112      CALL tra_unswap
113#endif
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#if defined key_obc || defined key_bdy || defined key_agrif
126      CALL tra_swap
127#endif
128 
129      ! set time step size (Euler/Leapfrog)
130      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
131      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
132      ENDIF
133
134      ! trends computation initialisation
135      IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter
136         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
137         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsn(:,:,:,jp_sal)
138      ENDIF
139
140      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
141         DO jn = 1, jpts
142            DO jk = 1, jpkm1
143               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
144            END DO
145         END DO
146      ELSE                                            ! Leap-Frog + Asselin filter time stepping
147         !
148         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
149         ELSE                 ;   CALL tra_nxt_fix( kt, 'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level
150         ENDIF
151      ENDIF 
152      !
153#if defined key_agrif
154      ! Update tracer at AGRIF zoom boundaries
155      CALL tra_unswap
156      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
157      CALL tra_swap
158#endif     
159      !
160      ! trends computation
161      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
162         DO jk = 1, jpkm1
163            zfact = 1.e0 / r2dtra(jk)             
164            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
165            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
166         END DO
167         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt )
168         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds )
169         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
170      END IF
171      !
172      !                        ! control print
173      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
174         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
175      !
176   END SUBROUTINE tra_nxt
177
178
179   SUBROUTINE tra_nxt_fix( kt, cdtype, ptb, ptn, pta, kjpt )
180      !!----------------------------------------------------------------------
181      !!                   ***  ROUTINE tra_nxt_fix  ***
182      !!
183      !! ** Purpose :   fixed volume: apply the Asselin time filter and
184      !!                swap the tracer fields.
185      !!
186      !! ** Method  : - Apply a Asselin time filter on now fields.
187      !!              - save in (ta,sa) an average over the three time levels
188      !!             which will be used to compute rdn and thus the semi-implicit
189      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
190      !!              - swap tracer fields to prepare the next time_step.
191      !!                This can be summurized for tempearture as:
192      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
193      !!             ztm = 0                                   otherwise
194      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
195      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
196      !!             tn  = ta 
197      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
198      !!
199      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
200      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
201      !!----------------------------------------------------------------------
202      INTEGER         , INTENT(in   )                               ::   kt       ! ocean 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 == nit000 )  THEN
215         IF(lwp) WRITE(numout,*)
216         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
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, cdtype, ptb, ptn, pta, 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      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
271      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
272      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
273      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
274      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
275      !!     
276      LOGICAL  ::   ll_tra, ll_tra_hpg, ll_traqsr   ! local logical
277      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
278      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
279      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
280      !!----------------------------------------------------------------------
281
282      IF( kt == nit000 ) THEN
283         IF(lwp) WRITE(numout,*)
284         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
285         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
286      ENDIF
287      !
288      IF( cdtype == 'TRA' )  THEN   
289         ll_tra     = .TRUE.           ! active tracers case 
290         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
291         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
292      ELSE                         
293         ll_tra     = .FALSE.          ! passive tracers case
294         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
295         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
296      ENDIF
297      !
298      DO jn = 1, kjpt     
299         DO jk = 1, jpkm1
300            zfact1 = atfp * rdttra(jk)
301            zfact2 = zfact1 / rau0
302            DO jj = 1, jpj
303               DO ji = 1, jpi
304                  ze3t_b = fse3t_b(ji,jj,jk)
305                  ze3t_n = fse3t_n(ji,jj,jk)
306                  ze3t_a = fse3t_a(ji,jj,jk)
307                  !                                         ! tracer content at Before, now and after
308                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
309                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
310                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
311                  !
312                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
313                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
314                  !
315                  ze3t_f = ze3t_n + atfp * ze3t_d
316                  ztc_f  = ztc_n  + atfp * ztc_d
317                  !
318                  IF( ll_tra .AND. jk == 1 ) THEN           ! first level only for T & S
319                      ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
320                      ztc_f  = ztc_f  - zfact1 * ( sbc_tsc(ji,jj,jn) - sbc_tsc_b(ji,jj,jn) )
321                  ENDIF
322                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
323                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
324
325                   ze3t_f = 1.e0 / ze3t_f
326                   ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
327                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
328                   !
329                   IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
330                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
331                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
332                   ENDIF
333               END DO
334            END DO
335         END DO
336         !
337      END DO
338      !
339   END SUBROUTINE tra_nxt_vvl
340
341   !!======================================================================
342END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.