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/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/TRP – NEMO

source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/TRP/trcatf.F90 @ 12443

Last change on this file since 12443 was 12443, checked in by davestorkey, 4 years ago

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