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/dev_r12377_KERNEL-06_techene_e3/src/TOP/TRP – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/TRP/trcatf.F90 @ 12779

Last change on this file since 12779 was 12779, checked in by techene, 4 years ago
modify TOP to be able to run with key_qco i.e. remove e3. gdep. and use a substitute instead (# include "domzgr_substitute.h90")
Property svn:keywords set to `Id`
File size: 17.2 KB

Rev	Line
[11057]	1	MODULE 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
[12779]	33	# if defined key_qco
	34	USE traatfqco
	35	# else
[11099]	36	USE traatf
[12779]	37	# endif
[7646]	38	USE bdy_oce , ONLY: ln_bdy
[6140]	39	USE trcbdy ! BDY open boundaries
[2528]	40	# if defined key_agrif
[941]	41	USE agrif_top_interp
[2528]	42	# endif
[9019]	43	!
	44	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	45	USE prtctl_trc ! Print control for debbuging
[941]	46
	47	IMPLICIT NONE
	48	PRIVATE
	49
[11057]	50	PUBLIC trc_atf ! routine called by step.F90
[2528]	51
[7881]	52	REAL(wp) :: rfact1, rfact2
[7872]	53
[12340]	54	!! * Substitutions
	55	# include "do_loop_substitute.h90"
[12779]	56	# include "domzgr_substitute.h90"
[941]	57	!!----------------------------------------------------------------------
[10067]	58	!! NEMO/TOP 4.0 , NEMO Consortium (2018)
[7753]	59	!! $Id$
[10068]	60	!! Software governed by the CeCILL license (see ./LICENSE)
[941]	61	!!----------------------------------------------------------------------
	62	CONTAINS
	63
[11057]	64	SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )
[941]	65	!!----------------------------------------------------------------------
[11057]	66	!! * ROUTINE trcatf *
[941]	67	!!
[11483]	68	!! ** Purpose : Apply the boundary condition on the after passive tracers fields and
	69	!! apply Asselin time filter to the now passive tracer fields if using leapfrog timestep
[941]	70	!!
[11057]	71	!! ** Method : Apply lateral boundary conditions on (uu(Kaa),vv(Kaa)) through
[941]	72	!! call to lbc_lnk routine
	73	!!
	74	!! For Arakawa or TVD Scheme :
[11483]	75	!! A Asselin time filter applied on now tracers tr(Kmm) to avoid
[941]	76	!! the divergence of two consecutive time-steps and tr arrays
	77	!! to prepare the next time_step:
[12724]	78	!! (tr(Kmm)) = (tr(Kmm)) + rn_atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
[941]	79	!!
	80	!!
[11483]	81	!! ** Action : - update tr(Kmm), tr(Kaa)
[941]	82	!!----------------------------------------------------------------------
[11057]	83	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	84	INTEGER , INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices
	85	REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers
[2715]	86	!
[7753]	87	INTEGER :: jk, jn ! dummy loop indices
[941]	88	REAL(wp) :: zfact ! temporary scalar
	89	CHARACTER (len=22) :: charout
[9019]	90	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztrdt ! 4D workspace
[941]	91	!!----------------------------------------------------------------------
[3294]	92	!
[11057]	93	IF( ln_timing ) CALL timing_start('trc_atf')
[3294]	94	!
	95	IF( kt == nittrc000 .AND. lwp ) THEN
[941]	96	WRITE(numout,*)
[11057]	97	WRITE(numout,*) 'trc_atf : Asselin time filtering on passive tracers'
[941]	98	ENDIF
[6140]	99	!
[5656]	100	#if defined key_agrif
	101	CALL Agrif_trc ! AGRIF zoom boundaries
	102	#endif
[9081]	103	! Update after tracer on domain lateral boundaries
[11057]	104	CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )
[941]	105
[11057]	106	IF( ln_bdy ) CALL trc_bdy( kt, Kbb, Kmm, Kaa )
[941]	107
[6140]	108	IF( l_trdtrc ) THEN ! trends: store now fields before the Asselin filter application
[9019]	109	ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
[10096]	110	ztrdt(:,:,:,:) = 0._wp
	111	IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend
	112	DO jn = 1, jptra
[11057]	113	CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
[12779]	114	END DO
[10096]	115	ENDIF
	116
	117	! total trend for the non-time-filtered variables.
[12724]	118	zfact = 1.0 / rn_Dt
[12779]	119	! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3ta*Ta)/e3tn; e3tn cancel from ts(Kmm) terms
[10096]	120	IF( ln_linssh ) THEN ! linear sea surface height only
	121	DO jn = 1, jptra
	122	DO jk = 1, jpkm1
[11057]	123	ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) zfact
[10096]	124	END DO
	125	END DO
	126	ELSE
	127	DO jn = 1, jptra
	128	DO jk = 1, jpkm1
[11057]	129	ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact
[10096]	130	END DO
	131	END DO
	132	ENDIF
	133	!
	134	DO jn = 1, jptra
[11057]	135	CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
[12779]	136	END DO
[10096]	137	!
	138	IF( ln_linssh ) THEN ! linear sea surface height only
	139	! Store now fields before applying the Asselin filter
	140	! in order to calculate Asselin filter trend later.
[11057]	141	ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm)
[10096]	142	ENDIF
	143
[2528]	144	ENDIF
[6140]	145	! ! Leap-Frog + Asselin filter time stepping
[12724]	146	IF( l_1st_euler .OR. ln_top_euler ) THEN ! Euler time-stepping
[11099]	147	!
	148	IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl
	149	! ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
	150	ztrdt(:,:,:,:) = 0._wp
	151	DO jn = 1, jptra
	152	CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
[12779]	153	END DO
[11099]	154	END IF
	155	!
	156	ELSE
[7881]	157	IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
[12779]	158	# if defined key_qco
	159	IF( ln_linssh ) THEN ; CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nittrc000, 'TRC', ptr, jptra ) ! linear ssh
	160	ELSE ; CALL tra_atf_qco_lf( kt, Kbb, Kmm, Kaa, nittrc000, rn_Dt, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
	161	# else
	162	IF( ln_linssh ) THEN ; CALL tra_atf_fix ( kt, Kbb, Kmm, Kaa, nittrc000, 'TRC', ptr, jptra ) ! linear ssh
	163	ELSE ; CALL tra_atf_vvl ( kt, Kbb, Kmm, Kaa, nittrc000, rn_Dt, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
	164	# endif
[7872]	165	ENDIF
	166	ELSE
[12779]	167	CALL trc_atf_off ( kt, Kbb, Kmm, Kaa, ptr ) ! offline
[2528]	168	ENDIF
[6140]	169	!
[11427]	170	CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
[2528]	171	ENDIF
[6140]	172	!
[10096]	173	IF( l_trdtrc .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt )
[2528]	174	DO jn = 1, jptra
	175	DO jk = 1, jpkm1
[12724]	176	zfact = 1._wp / rDt_trc
[11057]	177	ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact
[1175]	178	END DO
[11057]	179	CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
[2528]	180	END DO
	181	END IF
[10096]	182	IF( l_trdtrc ) DEALLOCATE( ztrdt )
[2528]	183	!
[12236]	184	IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
[941]	185	WRITE(charout, FMT="('nxt')")
	186	CALL prt_ctl_trc_info(charout)
[11057]	187	CALL prt_ctl_trc(tab4d=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
[941]	188	ENDIF
[2528]	189	!
[11057]	190	IF( ln_timing ) CALL timing_stop('trc_atf')
[3294]	191	!
[11057]	192	END SUBROUTINE trc_atf
[941]	193
[12779]	194	# if ! defined key_qco
[11057]	195	SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
[7872]	196	!!----------------------------------------------------------------------
[11057]	197	!! * ROUTINE tra_atf_off *
[7872]	198	!!
[11057]	199	!! !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
	200	!!
[7872]	201	!! ** Purpose : Time varying volume: apply the Asselin time filter
	202	!!
	203	!! ** Method : - Apply a thickness weighted Asselin time filter on now fields.
	204	!! - save in (ta,sa) a thickness weighted average over the three
	205	!! time levels which will be used to compute rdn and thus the semi-
	206	!! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
	207	!! - swap tracer fields to prepare the next time_step.
	208	!! This can be summurized for tempearture as:
	209	!! ztm = ( e3t_ntn + rbcp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] ) ln_dynhpg_imp = T
[12779]	210	!! /( e3t(:,:,jk,Kmm) + rbcp*[ e3t(:,:,jk,Kbb) - 2 e3t(:,:,jk,Kmm) + e3t(:,:,jk,Kaa) ] )
[7872]	211	!! ztm = 0 otherwise
[12724]	212	!! tb = ( e3t_ntn + rn_atfp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] )
[12779]	213	!! /( e3t(:,:,jk,Kmm) + rn_atfp*[ e3t(:,:,jk,Kbb) - 2 e3t(:,:,jk,Kmm) + e3t(:,:,jk,Kaa) ] )
[7872]	214	!! tn = ta
	215	!! ta = zt (NB: reset to 0 after eos_bn2 call)
	216	!!
	217	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	218	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	219	!!----------------------------------------------------------------------
[11057]	220	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	221	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices
	222	REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers
[7872]	223	!!
	224	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	225	REAL(wp) :: ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar
	226	REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d ! - -
	227	!!----------------------------------------------------------------------
	228	!
	229	IF( kt == nittrc000 ) THEN
	230	IF(lwp) WRITE(numout,*)
[11057]	231	IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
[7872]	232	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	233	IF( .NOT. ln_linssh ) THEN
[12724]	234	rfact1 = rn_atfp * rn_Dt
	235	rfact2 = rfact1 / rho0
[7872]	236	ENDIF
	237	!
	238	ENDIF
	239	!
	240	DO jn = 1, jptra
[12340]	241	DO_3D_11_11( 1, jpkm1 )
	242	ze3t_b = e3t(ji,jj,jk,Kbb)
	243	ze3t_n = e3t(ji,jj,jk,Kmm)
	244	ze3t_a = e3t(ji,jj,jk,Kaa)
	245	! ! tracer content at Before, now and after
	246	ztc_b = ptr(ji,jj,jk,jn,Kbb) * ze3t_b
	247	ztc_n = ptr(ji,jj,jk,jn,Kmm) * ze3t_n
	248	ztc_a = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
	249	!
	250	ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
	251	ztc_d = ztc_a - 2. * ztc_n + ztc_b
	252	!
[12724]	253	ze3t_f = ze3t_n + rn_atfp * ze3t_d
	254	ztc_f = ztc_n + rn_atfp * ztc_d
[12340]	255	!
	256	IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN ! first level
	257	ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
	258	ztc_f = ztc_f - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
	259	ENDIF
[7872]	260
[12340]	261	ze3t_f = 1.e0 / ze3t_f
	262	ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f ! time filtered "now" field
	263	!
	264	END_3D
[7872]	265	!
	266	END DO
	267	!
[11057]	268	END SUBROUTINE trc_atf_off
[12779]	269	# else
	270	SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
	271	!!----------------------------------------------------------------------
	272	!! * ROUTINE tra_atf_off *
	273	!!
	274	!! !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
	275	!!
	276	!! ** Purpose : Time varying volume: apply the Asselin time filter
	277	!!
	278	!! ** Method : - Apply a thickness weighted Asselin time filter on now fields.
	279	!! - save in (ta,sa) a thickness weighted average over the three
	280	!! time levels which will be used to compute rdn and thus the semi-
	281	!! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
	282	!! - swap tracer fields to prepare the next time_step.
	283	!! This can be summurized for tempearture as:
	284	!! ztm = ( e3t_ntn + rbcp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] ) ln_dynhpg_imp = T
	285	!! /( e3t(:,:,jk,Kmm) + rbcp*[ e3t(:,:,jk,Kbb) - 2 e3t(:,:,jk,Kmm) + e3t(:,:,jk,Kaa) ] )
	286	!! ztm = 0 otherwise
	287	!! tb = ( e3t_ntn + rn_atfp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] )
	288	!! /( e3t(:,:,jk,Kmm) + rn_atfp*[ e3t(:,:,jk,Kbb) - 2 e3t(:,:,jk,Kmm) + e3t(:,:,jk,Kaa) ] )
	289	!! tn = ta
	290	!! ta = zt (NB: reset to 0 after eos_bn2 call)
	291	!!
	292	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	293	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	294	!!----------------------------------------------------------------------
	295	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	296	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices
	297	REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers
	298	!!
	299	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	300	REAL(wp) :: ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar
	301	REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - -
	302	!!----------------------------------------------------------------------
	303	!
	304	IF( kt == nittrc000 ) THEN
	305	IF(lwp) WRITE(numout,*)
	306	IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
	307	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	308	IF( .NOT. ln_linssh ) THEN
	309	rfact1 = rn_atfp * rn_Dt
	310	rfact2 = rfact1 / rho0
	311	ENDIF
	312	!
	313	ENDIF
	314	!
	315	DO jn = 1, jptra
	316	DO_3D_11_11( 1, jpkm1 )
	317	ze3t_b = 1._wp + r3t(ji,jj,Kbb) * tmask(ji,jj,jk)
	318	ze3t_n = 1._wp + r3t(ji,jj,Kmm) * tmask(ji,jj,jk)
	319	ze3t_a = 1._wp + r3t(ji,jj,Kaa) * tmask(ji,jj,jk)
	320	! ! tracer content at Before, now and after
	321	ztc_b = ptr(ji,jj,jk,jn,Kbb) * ze3t_b
	322	ztc_n = ptr(ji,jj,jk,jn,Kmm) * ze3t_n
	323	ztc_a = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
	324	!
	325	ztc_d = ztc_a - 2. * ztc_n + ztc_b
	326	!
	327	ze3t_f = 1._wp + r3t_f(ji,jj)*tmask(ji,jj,jk)
	328	ztc_f = ztc_n + rn_atfp * ztc_d
	329	!
	330	IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN ! first level
	331	ztc_f = ztc_f - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
	332	ENDIF
[7872]	333
[12779]	334	ze3t_f = 1.e0 / ze3t_f
	335	ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f ! time filtered "now" field
	336	!
	337	END_3D
	338	!
	339	END DO
	340	!
	341	END SUBROUTINE trc_atf_off
	342	# endif
	343
[941]	344	#else
	345	!!----------------------------------------------------------------------
	346	!! Default option Empty module
	347	!!----------------------------------------------------------------------
[11099]	348	USE par_oce
	349	USE par_trc
[941]	350	CONTAINS
[11057]	351	SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )
[11099]	352	INTEGER , INTENT(in) :: kt
	353	INTEGER, INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices
	354	REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation
[11057]	355	WRITE(,) 'trc_atf: You should not have seen this print! error?', kt
	356	END SUBROUTINE trc_atf
[941]	357	#endif
	358	!!======================================================================
[11057]	359	END MODULE trcatf

Note: See TracBrowser for help on using the repository browser.

Download in other formats: