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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 5645

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

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

Download in other formats: