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.

tranxt.F90 in branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 5447

Last change on this file since 5447 was 5447, checked in by davestorkey, 9 years ago
Update UKMO/dev_r5107_mld_zint branch to revision 5442 of the trunk.
File size: 17.5 KB

Rev	Line
[3]	1	MODULE tranxt
	2	!!======================================================================
	3	!! * MODULE tranxt *
	4	!! Ocean active tracers: time stepping on temperature and salinity
	5	!!======================================================================
[1110]	6	!! History : OPA ! 1991-11 (G. Madec) Original code
	7	!! 7.0 ! 1993-03 (M. Guyon) symetrical conditions
	8	!! 8.0 ! 1996-02 (G. Madec & M. Imbard) opa release 8.0
	9	!! - ! 1996-04 (A. Weaver) Euler forward step
	10	!! 8.2 ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure grad.
	11	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	12	!! - ! 2002-11 (C. Talandier, A-M Treguier) Open boundaries
	13	!! - ! 2005-04 (C. Deltel) Add Asselin trend in the ML budget
	14	!! 2.0 ! 2006-02 (L. Debreu, C. Mazauric) Agrif implementation
	15	!! 3.0 ! 2008-06 (G. Madec) time stepping always done in trazdf
[1438]	16	!! 3.1 ! 2009-02 (G. Madec, R. Benshila) re-introduce the vvl option
[2528]	17	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) semi-implicit hpg with asselin filter + modified LF-RA
	18	!! - ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
[3]	19	!!----------------------------------------------------------------------
[503]	20
	21	!!----------------------------------------------------------------------
[2528]	22	!! tra_nxt : time stepping on tracers
	23	!! tra_nxt_fix : time stepping on tracers : fixed volume case
	24	!! tra_nxt_vvl : time stepping on tracers : variable volume case
[3]	25	!!----------------------------------------------------------------------
	26	USE oce ! ocean dynamics and tracers variables
	27	USE dom_oce ! ocean space and time domain variables
[2528]	28	USE sbc_oce ! surface boundary condition: ocean
[4990]	29	USE zdf_oce ! ocean vertical mixing
[1438]	30	USE domvvl ! variable volume
[1601]	31	USE dynspg_oce ! surface pressure gradient variables
	32	USE dynhpg ! hydrostatic pressure gradient
[4990]	33	USE trd_oce ! trends: ocean variables
	34	USE trdtra ! trends manager: tracers
	35	USE traqsr ! penetrative solar radiation (needed for nksr)
	36	USE phycst ! physical constant
	37	USE ldftra_oce ! lateral physics on tracers
	38	USE bdy_oce ! BDY open boundary condition variables
[3294]	39	USE bdytra ! open boundary condition (bdy_tra routine)
[4990]	40	!
[3]	41	USE in_out_manager ! I/O manager
	42	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[258]	43	USE prtctl ! Print control
[4990]	44	USE wrk_nemo ! Memory allocation
	45	USE timing ! Timing
[2528]	46	#if defined key_agrif
[389]	47	USE agrif_opa_update
	48	USE agrif_opa_interp
[2528]	49	#endif
[3]	50
	51	IMPLICIT NONE
	52	PRIVATE
	53
[2528]	54	PUBLIC tra_nxt ! routine called by step.F90
	55	PUBLIC tra_nxt_fix ! to be used in trcnxt
	56	PUBLIC tra_nxt_vvl ! to be used in trcnxt
[592]	57
[2715]	58	REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg
[1438]	59
[592]	60	!! * Substitutions
	61	# include "domzgr_substitute.h90"
[3]	62	!!----------------------------------------------------------------------
[2528]	63	!! NEMO/OPA 3.3 , NEMO-Consortium (2010)
[2715]	64	!! $Id$
	65	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	66	!!----------------------------------------------------------------------
	67	CONTAINS
	68
	69	SUBROUTINE tra_nxt( kt )
	70	!!----------------------------------------------------------------------
	71	!! * ROUTINE tranxt *
	72	!!
[1110]	73	!! ** Purpose : Apply the boundary condition on the after temperature
	74	!! and salinity fields, achieved the time stepping by adding
	75	!! the Asselin filter on now fields and swapping the fields.
[3]	76	!!
[1110]	77	!! ** Method : At this stage of the computation, ta and sa are the
	78	!! after temperature and salinity as the time stepping has
	79	!! been performed in trazdf_imp or trazdf_exp module.
[3]	80	!!
[1110]	81	!! - Apply lateral boundary conditions on (ta,sa)
	82	!! at the local domain boundaries through lbc_lnk call,
[4328]	83	!! at the one-way open boundaries (lk_bdy=T),
[4990]	84	!! at the AGRIF zoom boundaries (lk_agrif=T)
[1110]	85	!!
[1438]	86	!! - Update lateral boundary conditions on AGRIF children
	87	!! domains (lk_agrif=T)
[1110]	88	!!
	89	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	90	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
[503]	91	!!----------------------------------------------------------------------
	92	INTEGER, INTENT(in) :: kt ! ocean time-step index
	93	!!
[2528]	94	INTEGER :: jk, jn ! dummy loop indices
	95	REAL(wp) :: zfact ! local scalars
[3294]	96	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds
[3]	97	!!----------------------------------------------------------------------
[3294]	98	!
	99	IF( nn_timing == 1 ) CALL timing_start( 'tra_nxt')
	100	!
[1110]	101	IF( kt == nit000 ) THEN
	102	IF(lwp) WRITE(numout,*)
	103	IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
	104	IF(lwp) WRITE(numout,*) '~~~~~~~'
[2528]	105	!
[4230]	106	rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp) ! Brown & Campana parameter for semi-implicit hpg
[592]	107	ENDIF
	108
[1110]	109	! Update after tracer on domain lateral boundaries
	110	!
[4230]	111	CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp ) ! local domain boundaries (T-point, unchanged sign)
	112	CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
[1110]	113	!
[2528]	114	#if defined key_bdy
[3294]	115	IF( lk_bdy ) CALL bdy_tra( kt ) ! BDY open boundaries
[1110]	116	#endif
[392]	117	#if defined key_agrif
[2528]	118	CALL Agrif_tra ! AGRIF zoom boundaries
[389]	119	#endif
[1438]	120
	121	! set time step size (Euler/Leapfrog)
[2715]	122	IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dtra(:) = rdttra(:) ! at nit000 (Euler)
[4230]	123	ELSEIF( kt <= nit000 + 1 ) THEN ; r2dtra(:) = 2._wp* rdttra(:) ! at nit000 or nit000+1 (Leapfrog)
[1438]	124	ENDIF
[3]	125
[1110]	126	! trends computation initialisation
[2528]	127	IF( l_trdtra ) THEN ! store now fields before applying the Asselin filter
[3294]	128	CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
	129	ztrdt(:,:,:) = tsn(:,:,:,jp_tem)
	130	ztrds(:,:,:) = tsn(:,:,:,jp_sal)
[4990]	131	IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend
	132	CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
	133	CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
	134	ENDIF
[1110]	135	ENDIF
	136
[2528]	137	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step (only swap)
	138	DO jn = 1, jpts
	139	DO jk = 1, jpkm1
	140	tsn(:,:,jk,jn) = tsa(:,:,jk,jn)
	141	END DO
	142	END DO
	143	ELSE ! Leap-Frog + Asselin filter time stepping
	144	!
[5447]	145	IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa, &
	146	& sbc_tsc, sbc_tsc_b, jpts ) ! variable volume level (vvl)
	147	ELSE ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! fixed volume level
[2528]	148	ENDIF
	149	ENDIF
[2715]	150	!
[1438]	151	#if defined key_agrif
	152	! Update tracer at AGRIF zoom boundaries
	153	IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! children only
	154	#endif
[2715]	155	!
[1438]	156	! trends computation
[2528]	157	IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt
[1110]	158	DO jk = 1, jpkm1
[4990]	159	zfact = 1._wp / r2dtra(jk)
[2528]	160	ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
	161	ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
[1110]	162	END DO
[4990]	163	CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
	164	CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
[3294]	165	CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
[1438]	166	END IF
[2715]	167	!
[1438]	168	! ! control print
[2528]	169	IF(ln_ctl) CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt - Tn: ', mask1=tmask, &
	170	& tab3d_2=tsn(:,:,:,jp_sal), clinfo2= ' Sn: ', mask2=tmask )
[1438]	171	!
[4990]	172	IF( nn_timing == 1 ) CALL timing_stop('tra_nxt')
[3294]	173	!
[1438]	174	END SUBROUTINE tra_nxt
	175
	176
[3294]	177	SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
[1438]	178	!!----------------------------------------------------------------------
	179	!! * ROUTINE tra_nxt_fix *
	180	!!
	181	!! ** Purpose : fixed volume: apply the Asselin time filter and
	182	!! swap the tracer fields.
	183	!!
	184	!! ** Method : - Apply a Asselin time filter on now fields.
	185	!! - save in (ta,sa) an average over the three time levels
	186	!! which will be used to compute rdn and thus the semi-implicit
	187	!! hydrostatic pressure gradient (ln_dynhpg_imp = T)
	188	!! - swap tracer fields to prepare the next time_step.
	189	!! This can be summurized for tempearture as:
[2528]	190	!! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T
	191	!! ztm = 0 otherwise
	192	!! with rbcp=1/4 * (1-atfp^4) / (1-atfp)
[1438]	193	!! tb = tn + atfp*[ tb - 2 tn + ta ]
[2528]	194	!! tn = ta
[1438]	195	!! ta = ztm (NB: reset to 0 after eos_bn2 call)
	196	!!
	197	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	198	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	199	!!----------------------------------------------------------------------
[2715]	200	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	201	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2715]	202	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	203	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	204	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields
	205	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields
	206	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
	207	!
[2528]	208	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	209	LOGICAL :: ll_tra_hpg ! local logical
	210	REAL(wp) :: ztn, ztd ! local scalars
[1438]	211	!!----------------------------------------------------------------------
	212
[3294]	213	IF( kt == kit000 ) THEN
[1438]	214	IF(lwp) WRITE(numout,*)
[3294]	215	IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
[1438]	216	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	217	ENDIF
	218	!
[2528]	219	IF( cdtype == 'TRA' ) THEN ; ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg
	220	ELSE ; ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg
	221	ENDIF
	222	!
	223	DO jn = 1, kjpt
[1438]	224	!
[2528]	225	DO jk = 1, jpkm1
	226	DO jj = 1, jpj
	227	DO ji = 1, jpi
	228	ztn = ptn(ji,jj,jk,jn)
	229	ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers
	230	!
	231	ptb(ji,jj,jk,jn) = ztn + atfp * ztd ! ptb <-- filtered ptn
	232	ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta
	233	!
	234	IF( ll_tra_hpg ) pta(ji,jj,jk,jn) = ztn + rbcp * ztd ! pta <-- Brown & Campana average
[3]	235	END DO
[2528]	236	END DO
	237	END DO
[1110]	238	!
[2528]	239	END DO
[1438]	240	!
	241	END SUBROUTINE tra_nxt_fix
[3]	242
[1110]	243
[5447]	244	SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
[1438]	245	!!----------------------------------------------------------------------
	246	!! * ROUTINE tra_nxt_vvl *
	247	!!
	248	!! ** Purpose : Time varying volume: apply the Asselin time filter
	249	!! and swap the tracer fields.
	250	!!
	251	!! ** Method : - Apply a thickness weighted Asselin time filter on now fields.
	252	!! - save in (ta,sa) a thickness weighted average over the three
	253	!! time levels which will be used to compute rdn and thus the semi-
	254	!! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
	255	!! - swap tracer fields to prepare the next time_step.
	256	!! This can be summurized for tempearture as:
[2528]	257	!! ztm = ( e3t_ntn + rbcp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] ) ln_dynhpg_imp = T
	258	!! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] )
	259	!! ztm = 0 otherwise
[1438]	260	!! tb = ( e3t_ntn + atfp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] )
	261	!! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] )
	262	!! tn = ta
	263	!! ta = zt (NB: reset to 0 after eos_bn2 call)
	264	!!
	265	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	266	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	267	!!----------------------------------------------------------------------
[5447]	268	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	269	INTEGER , INTENT(in ) :: kit000 ! first time step index
	270	REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! time-step
	271	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	272	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	273	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields
	274	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields
	275	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
	276	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc ! surface tracer content
	277	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc_b ! before surface tracer content
	278
[1438]	279	!!
[5447]	280	LOGICAL :: ll_tra_hpg, ll_traqsr ! local logical
[2528]	281	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	282	REAL(wp) :: zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar
	283	REAL(wp) :: zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d ! - -
[1438]	284	!!----------------------------------------------------------------------
[3294]	285	!
	286	IF( kt == kit000 ) THEN
[1438]	287	IF(lwp) WRITE(numout,*)
[3294]	288	IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
[1438]	289	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	290	ENDIF
[2528]	291	!
	292	IF( cdtype == 'TRA' ) THEN
	293	ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg
	294	ll_traqsr = ln_traqsr ! active tracers case and solar penetration
	295	ELSE
	296	ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg
	297	ll_traqsr = .FALSE. ! active tracers case and NO solar penetration
	298	ENDIF
	299	!
	300	DO jn = 1, kjpt
	301	DO jk = 1, jpkm1
[5447]	302	zfact1 = atfp * p2dt(jk)
[2528]	303	zfact2 = zfact1 / rau0
	304	DO jj = 1, jpj
	305	DO ji = 1, jpi
	306	ze3t_b = fse3t_b(ji,jj,jk)
	307	ze3t_n = fse3t_n(ji,jj,jk)
	308	ze3t_a = fse3t_a(ji,jj,jk)
	309	! ! tracer content at Before, now and after
	310	ztc_b = ptb(ji,jj,jk,jn) * ze3t_b
	311	ztc_n = ptn(ji,jj,jk,jn) * ze3t_n
	312	ztc_a = pta(ji,jj,jk,jn) * ze3t_a
	313	!
	314	ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
	315	ztc_d = ztc_a - 2. * ztc_n + ztc_b
	316	!
	317	ze3t_f = ze3t_n + atfp * ze3t_d
	318	ztc_f = ztc_n + atfp * ztc_d
	319	!
[5447]	320	IF( jk == 1 ) THEN ! first level
	321	ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
	322	ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
[2528]	323	ENDIF
	324	IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & ! solar penetration (temperature only)
	325	& ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )
[1438]	326
[2528]	327	ze3t_f = 1.e0 / ze3t_f
	328	ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered
	329	ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta
	330	!
	331	IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only)
	332	ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d )
	333	pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average
	334	ENDIF
[1438]	335	END DO
	336	END DO
[2528]	337	END DO
	338	!
	339	END DO
[503]	340	!
[1438]	341	END SUBROUTINE tra_nxt_vvl
[3]	342
	343	!!======================================================================
	344	END MODULE tranxt

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

Download in other formats: