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.
trcnxt.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcnxt.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • Property svn:keywords set to Id
File size: 12.4 KB
Line 
1MODULE trcnxt
2   !!======================================================================
3   !!                       ***  MODULE  trcnxt  ***
4   !! Ocean passive tracers:  time stepping on passives tracers
5   !!======================================================================
6   !! History :  7.0  !  1991-11  (G. Madec)  Original code
7   !!                 !  1993-03  (M. Guyon)  symetrical conditions
8   !!                 !  1995-02  (M. Levy)   passive tracers
9   !!                 !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
10   !!            8.0  !  1996-04  (A. Weaver)  Euler forward step
11   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
12   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
13   !!                 !  2002-08  (G. Madec)  F90: Free form and module
14   !!                 !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
15   !!                 !  2004-03  (C. Ethe) passive tracers
16   !!                 !  2007-02  (C. Deltel) Diagnose ML trends for passive tracers
17   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
18   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
19   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
20   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
21   !!----------------------------------------------------------------------
22#if defined key_top
23   !!----------------------------------------------------------------------
24   !!   'key_top'                                                TOP models
25   !!----------------------------------------------------------------------
26   !!   trc_nxt     : time stepping on passive tracers
27   !!----------------------------------------------------------------------
28   USE oce_trc         ! ocean dynamics and tracers variables
29   USE trc             ! ocean passive tracers variables
30   USE trd_oce
31   USE trdtra
32   USE tranxt
33   USE bdy_oce   , ONLY: ln_bdy
34   USE trcbdy          ! BDY open boundaries
35# if defined key_agrif
36   USE agrif_top_interp
37# endif
38   !
39   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
40   USE prtctl_trc      ! Print control for debbuging
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   trc_nxt   ! routine called by step.F90
46
47   REAL(wp) ::   rfact1, rfact2
48
49   !!----------------------------------------------------------------------
50   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE trc_nxt( kt, Kmm, Krhs )
57      !!----------------------------------------------------------------------
58      !!                   ***  ROUTINE trcnxt  ***
59      !!
60      !! ** Purpose :   Compute the passive tracers fields at the
61      !!      next time-step from their temporal trends and swap the fields.
62      !!
63      !! ** Method  :   Apply lateral boundary conditions on (ua,va) through
64      !!      call to lbc_lnk routine
65      !!   default:
66      !!      arrays swap
67      !!         (trn) = (tra) ; (tra) = (0,0)
68      !!         (trb) = (trn)
69      !!
70      !!   For Arakawa or TVD Scheme :
71      !!      A Asselin time filter applied on now tracers (trn) to avoid
72      !!      the divergence of two consecutive time-steps and tr arrays
73      !!      to prepare the next time_step:
74      !!         (trb) = (trn) + atfp [ (trb) + (tra) - 2 (trn) ]
75      !!         (trn) = (tra) ; (tra) = (0,0)
76      !!
77      !!
78      !! ** Action  : - update trb, trn
79      !!----------------------------------------------------------------------
80      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
81      INTEGER, INTENT( in ) ::   Kmm, Krhs  ! time level indices
82      !
83      INTEGER  ::   jk, jn   ! dummy loop indices
84      REAL(wp) ::   zfact            ! temporary scalar
85      CHARACTER (len=22) :: charout
86      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace
87      !!----------------------------------------------------------------------
88      !
89      IF( ln_timing )   CALL timing_start('trc_nxt')
90      !
91      IF( kt == nittrc000 .AND. lwp ) THEN
92         WRITE(numout,*)
93         WRITE(numout,*) 'trc_nxt : time stepping on passive tracers'
94      ENDIF
95      !
96#if defined key_agrif
97      CALL Agrif_trc                   ! AGRIF zoom boundaries
98#endif
99      ! Update after tracer on domain lateral boundaries
100      CALL lbc_lnk( 'trcnxt', tra(:,:,:,:), 'T', 1. )   
101
102      IF( ln_bdy )  CALL trc_bdy( kt )
103
104      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
105         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
106         ztrdt(:,:,:,:)  = 0._wp
107         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend
108            DO jn = 1, jptra
109               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
110            ENDDO
111         ENDIF
112
113         ! total trend for the non-time-filtered variables.
114         zfact = 1.0 / rdttrc
115         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms
116         IF( ln_linssh ) THEN       ! linear sea surface height only
117            DO jn = 1, jptra
118               DO jk = 1, jpkm1
119                  ztrdt(:,:,jk,jn) = ( tra(:,:,jk,jn)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - trn(:,:,jk,jn)) * zfact
120               END DO
121            END DO
122         ELSE
123            DO jn = 1, jptra
124               DO jk = 1, jpkm1
125                  ztrdt(:,:,jk,jn) = ( tra(:,:,jk,jn) - trn(:,:,jk,jn) ) * zfact
126               END DO
127            END DO
128         ENDIF
129         !
130         DO jn = 1, jptra
131            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
132         ENDDO
133         !
134         IF( ln_linssh ) THEN       ! linear sea surface height only
135            ! Store now fields before applying the Asselin filter
136            ! in order to calculate Asselin filter trend later.
137            ztrdt(:,:,:,:) = trn(:,:,:,:) 
138         ENDIF
139
140      ENDIF
141      !                                ! Leap-Frog + Asselin filter time stepping
142      IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping (only swap)
143         DO jn = 1, jptra
144            DO jk = 1, jpkm1
145               trn(:,:,jk,jn) = tra(:,:,jk,jn)
146               trb(:,:,jk,jn) = trn(:,:,jk,jn) 
147            END DO
148         END DO
149         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
150            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
151            ztrdt(:,:,:,:) = 0._wp           
152            DO jn = 1, jptra
153               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
154            ENDDO
155         END IF
156         !
157      ELSE     
158         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
159            IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nittrc000,         'TRC', trb, trn, tra, jptra )  !     linear ssh
160            ELSE                   ;   CALL tra_nxt_vvl( kt, Kmm, Krhs, nittrc000, rdttrc, 'TRC', trb, trn, tra,      &
161              &                                                                   sbc_trc, sbc_trc_b, jptra )  ! non-linear ssh
162            ENDIF
163         ELSE
164                                       CALL trc_nxt_off( kt )       ! offline
165         ENDIF
166         !
167         CALL lbc_lnk_multi( 'trcnxt', trb(:,:,:,:), 'T', 1._wp, trn(:,:,:,:), 'T', 1._wp, tra(:,:,:,:), 'T', 1._wp )
168      ENDIF
169      !
170      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
171         DO jn = 1, jptra
172            DO jk = 1, jpkm1
173               zfact = 1._wp / r2dttrc 
174               ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact 
175            END DO
176            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
177         END DO
178      END IF
179      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
180      !
181      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
182         WRITE(charout, FMT="('nxt')")
183         CALL prt_ctl_trc_info(charout)
184         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
185      ENDIF
186      !
187      IF( ln_timing )   CALL timing_stop('trc_nxt')
188      !
189   END SUBROUTINE trc_nxt
190
191
192   SUBROUTINE trc_nxt_off( kt )
193      !!----------------------------------------------------------------------
194      !!                   ***  ROUTINE tra_nxt_vvl  ***
195      !!
196      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
197      !!                and swap the tracer fields.
198      !!
199      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
200      !!              - save in (ta,sa) a thickness weighted average over the three
201      !!             time levels which will be used to compute rdn and thus the semi-
202      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
203      !!              - swap tracer fields to prepare the next time_step.
204      !!                This can be summurized for tempearture as:
205      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
206      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
207      !!             ztm = 0                                                       otherwise
208      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
209      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
210      !!             tn  = ta
211      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
212      !!
213      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
214      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
215      !!----------------------------------------------------------------------
216      INTEGER , INTENT(in   )   ::  kt       ! ocean time-step index
217      !!     
218      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
219      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
220      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
221      !!----------------------------------------------------------------------
222      !
223      IF( kt == nittrc000 )  THEN
224         IF(lwp) WRITE(numout,*)
225         IF(lwp) WRITE(numout,*) 'trc_nxt_off : time stepping'
226         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
227         IF( .NOT. ln_linssh ) THEN
228            rfact1 = atfp * rdttrc
229            rfact2 = rfact1 / rau0
230         ENDIF
231        ! 
232      ENDIF
233      !
234      DO jn = 1, jptra     
235         DO jk = 1, jpkm1
236            DO jj = 1, jpj
237               DO ji = 1, jpi
238                  ze3t_b = e3t_b(ji,jj,jk)
239                  ze3t_n = e3t_n(ji,jj,jk)
240                  ze3t_a = e3t_a(ji,jj,jk)
241                  !                                         ! tracer content at Before, now and after
242                  ztc_b  = trb(ji,jj,jk,jn) * ze3t_b
243                  ztc_n  = trn(ji,jj,jk,jn) * ze3t_n
244                  ztc_a  = tra(ji,jj,jk,jn) * ze3t_a
245                  !
246                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
247                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
248                  !
249                  ze3t_f = ze3t_n + atfp * ze3t_d
250                  ztc_f  = ztc_n  + atfp * ztc_d
251                  !
252                  IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
253                     ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
254                     ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
255                  ENDIF
256
257                  ze3t_f = 1.e0 / ze3t_f
258                  trb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
259                  trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn)     ! ptn <-- pta
260                  !
261               END DO
262            END DO
263         END DO
264         !
265      END DO
266      !
267   END SUBROUTINE trc_nxt_off
268
269#else
270   !!----------------------------------------------------------------------
271   !!   Default option                                         Empty module
272   !!----------------------------------------------------------------------
273CONTAINS
274   SUBROUTINE trc_nxt( kt ) 
275      INTEGER, INTENT(in) :: kt
276      WRITE(*,*) 'trc_nxt: You should not have seen this print! error?', kt
277   END SUBROUTINE trc_nxt
278#endif
279   !!======================================================================
280END MODULE trcnxt
Note: See TracBrowser for help on using the repository browser.