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.
trcatf.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/trcatf.F90 @ 11483

Last change on this file since 11483 was 11483, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Update more header comments and change logs.
Also remove debugging output from trcstp.F90 and correct bug in step_c1d.F90. Note C1D untested in this branch.
  • Property svn:keywords set to Id
File size: 13.1 KB
Line 
1MODULE trcatf
2   !!======================================================================
3   !!                       ***  MODULE  trcatf  ***
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   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename trcnxt.F90 -> trcatf.F90. Now only does time filtering.
22   !!----------------------------------------------------------------------
23#if defined key_top
24   !!----------------------------------------------------------------------
25   !!   'key_top'                                                TOP models
26   !!----------------------------------------------------------------------
27   !!   trc_atf     : time stepping on passive tracers
28   !!----------------------------------------------------------------------
29   USE oce_trc         ! ocean dynamics and tracers variables
30   USE trc             ! ocean passive tracers variables
31   USE trd_oce
32   USE trdtra
33   USE traatf
34   USE bdy_oce   , ONLY: ln_bdy
35   USE trcbdy          ! BDY open boundaries
36# if defined key_agrif
37   USE agrif_top_interp
38# endif
39   !
40   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
41   USE prtctl_trc      ! Print control for debbuging
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   trc_atf   ! routine called by step.F90
47
48   REAL(wp) ::   rfact1, rfact2
49
50   !!----------------------------------------------------------------------
51   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )
58      !!----------------------------------------------------------------------
59      !!                   ***  ROUTINE trcatf  ***
60      !!
61      !! ** Purpose :   Apply the boundary condition on the after passive tracers fields and
62      !!      apply Asselin time filter to the now passive tracer fields if using leapfrog timestep
63      !!
64      !! ** Method  :   Apply lateral boundary conditions on (uu(Kaa),vv(Kaa)) through
65      !!      call to lbc_lnk routine
66      !!
67      !!   For Arakawa or TVD Scheme :
68      !!      A Asselin time filter applied on now tracers tr(Kmm) to avoid
69      !!      the divergence of two consecutive time-steps and tr arrays
70      !!      to prepare the next time_step:
71      !!         (tr(Kmm)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
72      !!
73      !!
74      !! ** Action  : - update tr(Kmm), tr(Kaa)
75      !!----------------------------------------------------------------------
76      INTEGER                                   , INTENT( in )  :: kt             ! ocean time-step index
77      INTEGER                                   , INTENT( in )  :: Kbb, Kmm, Kaa ! time level indices
78      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers
79      !
80      INTEGER  ::   jk, jn   ! dummy loop indices
81      REAL(wp) ::   zfact            ! temporary scalar
82      CHARACTER (len=22) :: charout
83      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace
84      !!----------------------------------------------------------------------
85      !
86      IF( ln_timing )   CALL timing_start('trc_atf')
87      !
88      IF( kt == nittrc000 .AND. lwp ) THEN
89         WRITE(numout,*)
90         WRITE(numout,*) 'trc_atf : Asselin time filtering on passive tracers'
91      ENDIF
92      !
93#if defined key_agrif
94      CALL Agrif_trc                   ! AGRIF zoom boundaries
95#endif
96      ! Update after tracer on domain lateral boundaries
97      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )   
98
99      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Kaa )
100
101      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
102         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
103         ztrdt(:,:,:,:)  = 0._wp
104         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend
105            DO jn = 1, jptra
106               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
107            ENDDO
108         ENDIF
109
110         ! total trend for the non-time-filtered variables.
111         zfact = 1.0 / rdttrc
112         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms
113         IF( ln_linssh ) THEN       ! linear sea surface height only
114            DO jn = 1, jptra
115               DO jk = 1, jpkm1
116                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) * zfact
117               END DO
118            END DO
119         ELSE
120            DO jn = 1, jptra
121               DO jk = 1, jpkm1
122                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact
123               END DO
124            END DO
125         ENDIF
126         !
127         DO jn = 1, jptra
128            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
129         ENDDO
130         !
131         IF( ln_linssh ) THEN       ! linear sea surface height only
132            ! Store now fields before applying the Asselin filter
133            ! in order to calculate Asselin filter trend later.
134            ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm) 
135         ENDIF
136
137      ENDIF
138      !                                ! Leap-Frog + Asselin filter time stepping
139      IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping
140         !
141         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
142            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
143            ztrdt(:,:,:,:) = 0._wp           
144            DO jn = 1, jptra
145               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
146            ENDDO
147         END IF
148         !
149      ELSE     
150         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
151            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
152            ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rdttrc, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
153            ENDIF
154         ELSE
155                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )       ! offline
156         ENDIF
157         !
158         CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
159      ENDIF
160      !
161      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
162         DO jn = 1, jptra
163            DO jk = 1, jpkm1
164               zfact = 1._wp / r2dttrc 
165               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
166            END DO
167            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
168         END DO
169      END IF
170      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
171      !
172      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
173         WRITE(charout, FMT="('nxt')")
174         CALL prt_ctl_trc_info(charout)
175         CALL prt_ctl_trc(tab4d=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
176      ENDIF
177      !
178      IF( ln_timing )   CALL timing_stop('trc_atf')
179      !
180   END SUBROUTINE trc_atf
181
182
183   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
184      !!----------------------------------------------------------------------
185      !!                   ***  ROUTINE tra_atf_off  ***
186      !!
187      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
188      !!
189      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
190      !!
191      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
192      !!              - save in (ta,sa) a thickness weighted average over the three
193      !!             time levels which will be used to compute rdn and thus the semi-
194      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
195      !!              - swap tracer fields to prepare the next time_step.
196      !!                This can be summurized for tempearture as:
197      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
198      !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )   
199      !!             ztm = 0                                                       otherwise
200      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
201      !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )
202      !!             tn  = ta
203      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
204      !!
205      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
206      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
207      !!----------------------------------------------------------------------
208      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
209      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
210      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
211      !!     
212      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
213      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
214      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
215      !!----------------------------------------------------------------------
216      !
217      IF( kt == nittrc000 )  THEN
218         IF(lwp) WRITE(numout,*)
219         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
220         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
221         IF( .NOT. ln_linssh ) THEN
222            rfact1 = atfp * rdttrc
223            rfact2 = rfact1 / rau0
224         ENDIF
225        ! 
226      ENDIF
227      !
228      DO jn = 1, jptra     
229         DO jk = 1, jpkm1
230            DO jj = 1, jpj
231               DO ji = 1, jpi
232                  ze3t_b = e3t(ji,jj,jk,Kbb)
233                  ze3t_n = e3t(ji,jj,jk,Kmm)
234                  ze3t_a = e3t(ji,jj,jk,Kaa)
235                  !                                         ! tracer content at Before, now and after
236                  ztc_b  = ptr(ji,jj,jk,jn,Kbb)  * ze3t_b
237                  ztc_n  = ptr(ji,jj,jk,jn,Kmm)  * ze3t_n
238                  ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
239                  !
240                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
241                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
242                  !
243                  ze3t_f = ze3t_n + atfp * ze3t_d
244                  ztc_f  = ztc_n  + atfp * ztc_d
245                  !
246                  IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
247                     ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
248                     ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
249                  ENDIF
250
251                  ze3t_f = 1.e0 / ze3t_f
252                  ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
253                  !
254               END DO
255            END DO
256         END DO
257         !
258      END DO
259      !
260   END SUBROUTINE trc_atf_off
261
262#else
263   !!----------------------------------------------------------------------
264   !!   Default option                                         Empty module
265   !!----------------------------------------------------------------------
266   USE par_oce
267   USE par_trc
268CONTAINS
269   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr ) 
270      INTEGER                                   , INTENT(in)    :: kt
271      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm, Kaa ! time level indices
272      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr           ! passive tracers and RHS of tracer equation
273      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
274   END SUBROUTINE trc_atf
275#endif
276   !!======================================================================
277END MODULE trcatf
Note: See TracBrowser for help on using the repository browser.