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/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 5600

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

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

Download in other formats: