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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago
set proper svn properties to all files...
Property svn:keywords set to `Id`
File size: 17.3 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
[2148]	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	!!----------------------------------------------------------------------
[2148]	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
[2148]	28	USE sbc_oce ! surface boundary condition: ocean
[3]	29	USE zdf_oce ! ???
[1438]	30	USE domvvl ! variable volume
[1601]	31	USE dynspg_oce ! surface pressure gradient variables
	32	USE dynhpg ! hydrostatic pressure gradient
[2024]	33	USE trdmod_oce ! ocean space and time domain variables
	34	USE trdtra ! ocean active tracers trends
[888]	35	USE phycst
[2120]	36	USE obc_oce
[888]	37	USE obctra ! open boundary condition (obc_tra routine)
[911]	38	USE bdytra ! Unstructured open boundary condition (bdy_tra routine)
[3]	39	USE in_out_manager ! I/O manager
	40	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[258]	41	USE prtctl ! Print control
[2148]	42	USE traqsr ! penetrative solar radiation (needed for nksr)
[2024]	43	USE traswp ! swap array
[2236]	44	USE obc_oce
	45	#if defined key_agrif
[389]	46	USE agrif_opa_update
	47	USE agrif_opa_interp
[2236]	48	#endif
[3]	49
	50	IMPLICIT NONE
	51	PRIVATE
	52
[2024]	53	PUBLIC tra_nxt ! routine called by step.F90
	54	PUBLIC tra_nxt_fix ! to be used in trcnxt
	55	PUBLIC tra_nxt_vvl ! to be used in trcnxt
[592]	56
[2148]	57	REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg
[2082]	58	REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
[1438]	59
[592]	60	!! * Substitutions
	61	# include "domzgr_substitute.h90"
[3]	62	!!----------------------------------------------------------------------
[2272]	63	!! NEMO/OPA 3.3 , NEMO-Consortium (2010)
	64	!! $Id $
[2236]	65	!! Software governed by the CeCILL licence (NEMOGCM/License_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,
	83	!! at the radiative open boundaries (lk_obc=T),
	84	!! at the relaxed open boundaries (lk_bdy=T), and
	85	!! at the AGRIF zoom boundaries (lk_agrif=T)
	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	!!
[2257]	95	INTEGER :: jk, jn ! dummy loop indices
	96	REAL(wp) :: zfact ! local scalars
[2024]	97	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds
[3]	98	!!----------------------------------------------------------------------
	99
[1110]	100	IF( kt == nit000 ) THEN
	101	IF(lwp) WRITE(numout,*)
	102	IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
	103	IF(lwp) WRITE(numout,*) '~~~~~~~'
[2148]	104	!
	105	rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp) ! Brown & Campana parameter for semi-implicit hpg
[592]	106	ENDIF
	107
[1110]	108	! Update after tracer on domain lateral boundaries
	109	!
[2024]	110	CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ! local domain boundaries (T-point, unchanged sign)
	111	CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
[1110]	112	!
[2024]	113	#if defined key_obc \|\| defined key_bdy \|\| defined key_agrif
	114	CALL tra_unswap
	115	#endif
[2120]	116
[2148]	117	#if defined key_obc
[1970]	118	IF( lk_obc ) CALL obc_tra( kt ) ! OBC open boundaries
[3]	119	#endif
[2120]	120	#if defined key_bdy
	121	IF( lk_bdy ) CALL bdy_tra( kt ) ! BDY open boundaries
[1110]	122	#endif
[392]	123	#if defined key_agrif
[2272]	124	CALL Agrif_tra ! AGRIF zoom boundaries
[389]	125	#endif
[2120]	126
[2024]	127	#if defined key_obc \|\| defined key_bdy \|\| defined key_agrif
	128	CALL tra_swap
	129	#endif
[1438]	130
	131	! set time step size (Euler/Leapfrog)
[2082]	132	IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt(:) = rdttra(:) ! at nit000 (Euler)
	133	ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog)
[1438]	134	ENDIF
[3]	135
[1110]	136	! trends computation initialisation
[2272]	137	IF( l_trdtra ) THEN ! store now fields before applying the Asselin filter
[2024]	138	ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsn(:,:,:,jp_tem)
	139	ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsn(:,:,:,jp_sal)
[1110]	140	ENDIF
	141
[2272]	142	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step (only swap)
[2257]	143	DO jn = 1, jpts
	144	DO jk = 1, jpkm1
	145	tsn(:,:,jk,jn) = tsa(:,:,jk,jn)
	146	END DO
	147	END DO
[2272]	148	ELSE ! Leap-Frog + Asselin filter time stepping
	149	!
[2257]	150	IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts ) ! variable volume level (vvl)
	151	ELSE ; CALL tra_nxt_fix( kt, 'TRA', tsb, tsn, tsa, jpts ) ! fixed volume level
	152	ENDIF
	153	ENDIF
[1438]	154
	155	#if defined key_agrif
[2272]	156	! Update tracer at AGRIF zoom boundaries
[2024]	157	CALL tra_unswap
[1438]	158	IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! children only
[2024]	159	CALL tra_swap
[1438]	160	#endif
	161
	162	! trends computation
[2272]	163	IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt
[1110]	164	DO jk = 1, jpkm1
[2082]	165	zfact = 1.e0 / r2dt(jk)
[2024]	166	ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
	167	ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
[1110]	168	END DO
[2083]	169	CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt )
	170	CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds )
[2024]	171	DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds )
[1438]	172	END IF
	173
	174	! ! control print
[2024]	175	IF(ln_ctl) CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt - Tn: ', mask1=tmask, &
	176	& tab3d_2=tsn(:,:,:,jp_sal), clinfo2= ' Sn: ', mask2=tmask )
[1438]	177	!
	178	END SUBROUTINE tra_nxt
	179
[2104]	180
[2240]	181	SUBROUTINE tra_nxt_fix( kt, cdtype, ptb, ptn, pta, kjpt )
[1438]	182	!!----------------------------------------------------------------------
	183	!! * ROUTINE tra_nxt_fix *
	184	!!
	185	!! ** Purpose : fixed volume: apply the Asselin time filter and
	186	!! swap the tracer fields.
	187	!!
	188	!! ** Method : - Apply a Asselin time filter on now fields.
	189	!! - save in (ta,sa) an average over the three time levels
	190	!! which will be used to compute rdn and thus the semi-implicit
	191	!! hydrostatic pressure gradient (ln_dynhpg_imp = T)
	192	!! - swap tracer fields to prepare the next time_step.
	193	!! This can be summurized for tempearture as:
[2148]	194	!! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T
	195	!! ztm = 0 otherwise
	196	!! with rbcp=1/4 * (1-atfp^4) / (1-atfp)
[1438]	197	!! tb = tn + atfp*[ tb - 2 tn + ta ]
[2148]	198	!! tn = ta
[1438]	199	!! ta = ztm (NB: reset to 0 after eos_bn2 call)
	200	!!
	201	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	202	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	203	!!----------------------------------------------------------------------
[2240]	204	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	205	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	206	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	207	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields
	208	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields
	209	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
[2272]	210	!!
	211	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	212	LOGICAL :: ll_tra_hpg ! local logical
	213	REAL(wp) :: ztn, ztd, ztm ! local scalars
[1438]	214	!!----------------------------------------------------------------------
	215
[2104]	216	IF( kt == nit000 ) THEN
[1438]	217	IF(lwp) WRITE(numout,*)
	218	IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
	219	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	220	ENDIF
	221	!
[2272]	222	IF( cdtype == 'TRA' ) THEN ; ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg
	223	ELSE ; ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg
[2257]	224	ENDIF
[2120]	225	!
[2257]	226	DO jn = 1, kjpt
[2120]	227	!
[2257]	228	DO jk = 1, jpkm1
	229	DO jj = 1, jpj
	230	DO ji = 1, jpi
[2272]	231	IF( ll_tra_hpg ) ztn = ptn(ji,jj,jk,jn) ! semi-implicit hpg: keep tn, sn in memory
[2257]	232	!
	233	ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ! time laplacian on tracers
	234	!
	235	ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd ! ptb <-- filtered ptn
	236	ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta
	237	!
[2272]	238	IF( ll_tra_hpg ) pta(ji,jj,jk,jn) = ztn + rbcp * ztd ! pta <-- Brown & Campana average
[3]	239	END DO
[2257]	240	END DO
	241	END DO
[2024]	242	!
[2257]	243	END DO
[1438]	244	!
	245	END SUBROUTINE tra_nxt_fix
[3]	246
[2104]	247
[2148]	248	SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt )
[1438]	249	!!----------------------------------------------------------------------
	250	!! * ROUTINE tra_nxt_vvl *
	251	!!
	252	!! ** Purpose : Time varying volume: apply the Asselin time filter
	253	!! and swap the tracer fields.
	254	!!
	255	!! ** Method : - Apply a thickness weighted Asselin time filter on now fields.
	256	!! - save in (ta,sa) a thickness weighted average over the three
	257	!! time levels which will be used to compute rdn and thus the semi-
	258	!! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
	259	!! - swap tracer fields to prepare the next time_step.
	260	!! This can be summurized for tempearture as:
[2148]	261	!! ztm = ( e3t_ntn + rbcp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] ) ln_dynhpg_imp = T
	262	!! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] )
	263	!! ztm = 0 otherwise
[1438]	264	!! tb = ( e3t_ntn + atfp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] )
	265	!! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] )
	266	!! tn = ta
	267	!! ta = zt (NB: reset to 0 after eos_bn2 call)
	268	!!
	269	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	270	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	271	!!----------------------------------------------------------------------
[2148]	272	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	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
[1438]	278	!!
[2272]	279	LOGICAL :: ll_tra, ll_tra_hpg, ll_traqsr ! local logical
	280	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	281	REAL(wp) :: ztc_a , ztc_n , ztc_b ! local scalar
	282	REAL(wp) :: ztc_f , ztc_d , ztc_m ! - -
	283	REAL(wp) :: ze3t_b, ze3t_n, ze3t_a ! - -
	284	REAL(wp) :: ze3t_f, ze3t_d, ze3t_m ! - -
	285	REAL(wp) :: zfact1, zfact2 ! - -
[1438]	286	!!----------------------------------------------------------------------
	287
[2104]	288	IF( kt == nit000 ) THEN
[1438]	289	IF(lwp) WRITE(numout,*)
	290	IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
	291	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	292	ENDIF
[2024]	293	!
[2272]	294	IF( cdtype == 'TRA' ) THEN
	295	ll_tra = .TRUE. ! active tracers case
	296	ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg
	297	ll_traqsr = ln_traqsr ! active tracers case and semi-implicit hpg
	298	ELSE
	299	ll_tra = .FALSE. ! passive tracers case
	300	ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg
	301	ll_traqsr = .FALSE. ! active tracers case and semi-implicit hpg
[2257]	302	ENDIF
[2120]	303	!
[2257]	304	DO jn = 1, kjpt
	305	DO jk = 1, jpkm1
	306	zfact1 = atfp * rdttra(jk)
	307	zfact2 = zfact1 / rau0
	308	DO jj = 1, jpj
	309	DO ji = 1, jpi
	310	ze3t_b = fse3t_b(ji,jj,jk)
	311	ze3t_n = fse3t_n(ji,jj,jk)
	312	ze3t_a = fse3t_a(ji,jj,jk)
	313	! ! tracer content at Before, now and after
	314	ztc_b = ptb(ji,jj,jk,jn) * ze3t_b
	315	ztc_n = ptn(ji,jj,jk,jn) * ze3t_n
	316	ztc_a = pta(ji,jj,jk,jn) * ze3t_a
	317	!
	318	ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
	319	ztc_d = ztc_a - 2. * ztc_n + ztc_b
	320	!
	321	ze3t_f = ze3t_n + atfp * ze3t_d
	322	ztc_f = ztc_n + atfp * ztc_d
[2272]	323	!
	324	IF( ll_tra .AND. jk == 1 ) THEN ! first level only for T & S
[2257]	325	ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
	326	ztc_f = ztc_f - zfact1 * ( sbc_tsc(ji,jj,jn) - sbc_tsc_b(ji,jj,jn) )
	327	ENDIF
[2272]	328	IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & ! solar penetration (temperature only)
[2257]	329	& ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )
[2148]	330
[2257]	331	ze3t_f = 1.e0 / ze3t_f
	332	ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered
	333	ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta
	334	!
[2272]	335	IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only)
[2257]	336	ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d )
	337	pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average
	338	ENDIF
[1438]	339	END DO
[2257]	340	END DO
	341	END DO
	342	!
	343	END DO
[503]	344	!
[1438]	345	END SUBROUTINE tra_nxt_vvl
[3]	346
	347	!!======================================================================
	348	END MODULE tranxt

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

Download in other formats: