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_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcatf.F90 @ 13373

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

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

Download in other formats: