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

Context Navigation

source: branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90 @ 2091

Last change on this file since 2091 was 2091, checked in by mlelod, 14 years ago
ticket: #663 cosmetic changes
Property svn:eol-style set to `native` 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
[1870]	17	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) semi-implicit hpg with asselin filter + modified LF-RA
[3]	18	!!----------------------------------------------------------------------
[503]	19
	20	!!----------------------------------------------------------------------
[1110]	21	!! tra_nxt : time stepping on temperature and salinity
[1438]	22	!! tra_nxt_fix : time stepping on temperature and salinity : fixed volume case
	23	!! tra_nxt_vvl : time stepping on temperature and salinity : variable volume case
[3]	24	!!----------------------------------------------------------------------
	25	USE oce ! ocean dynamics and tracers variables
	26	USE dom_oce ! ocean space and time domain variables
[1975]	27	USE sbc_oce ! surface boundary condition: ocean
[3]	28	USE zdf_oce ! ???
[1438]	29	USE domvvl ! variable volume
[1601]	30	USE dynspg_oce ! surface pressure gradient variables
	31	USE dynhpg ! hydrostatic pressure gradient
[888]	32	USE trdmod_oce ! ocean variables trends
	33	USE trdmod ! ocean active tracers trends
	34	USE phycst
	35	USE obctra ! open boundary condition (obc_tra routine)
[911]	36	USE bdytra ! Unstructured open boundary condition (bdy_tra routine)
[3]	37	USE in_out_manager ! I/O manager
	38	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[258]	39	USE prtctl ! Print control
[1975]	40	USE traqsr ! penetrative solar radiation (needed for nksr)
[389]	41	USE agrif_opa_update
	42	USE agrif_opa_interp
[3]	43
	44	IMPLICIT NONE
	45	PRIVATE
	46
[1110]	47	PUBLIC tra_nxt ! routine called by step.F90
[592]	48
[1975]	49	REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg
[1847]	50	REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
[1438]	51
[592]	52	!! * Substitutions
	53	# include "domzgr_substitute.h90"
[3]	54	!!----------------------------------------------------------------------
[1847]	55	!! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
[1146]	56	!! $Id$
[503]	57	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	58	!!----------------------------------------------------------------------
	59
	60	CONTAINS
	61
	62	SUBROUTINE tra_nxt( kt )
	63	!!----------------------------------------------------------------------
	64	!! * ROUTINE tranxt *
	65	!!
[1110]	66	!! ** Purpose : Apply the boundary condition on the after temperature
	67	!! and salinity fields, achieved the time stepping by adding
	68	!! the Asselin filter on now fields and swapping the fields.
[3]	69	!!
[1110]	70	!! ** Method : At this stage of the computation, ta and sa are the
	71	!! after temperature and salinity as the time stepping has
	72	!! been performed in trazdf_imp or trazdf_exp module.
[3]	73	!!
[1110]	74	!! - Apply lateral boundary conditions on (ta,sa)
	75	!! at the local domain boundaries through lbc_lnk call,
	76	!! at the radiative open boundaries (lk_obc=T),
	77	!! at the relaxed open boundaries (lk_bdy=T), and
	78	!! at the AGRIF zoom boundaries (lk_agrif=T)
	79	!!
[1438]	80	!! - Update lateral boundary conditions on AGRIF children
	81	!! domains (lk_agrif=T)
[1110]	82	!!
	83	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	84	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
[1870]	85	!!
	86	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[503]	87	!!----------------------------------------------------------------------
	88	USE oce, ONLY : ztrdt => ua ! use ua as 3D workspace
	89	USE oce, ONLY : ztrds => va ! use va as 3D workspace
[3]	90	!!
[503]	91	INTEGER, INTENT(in) :: kt ! ocean time-step index
	92	!!
[1847]	93	INTEGER :: jk ! dummy loop indices
	94	REAL(wp) :: zfact ! temporary scalars
[3]	95	!!----------------------------------------------------------------------
	96
[1870]	97	IF( kt == nit000 ) THEN !== initialisation ==!
[1110]	98	IF(lwp) WRITE(numout,*)
	99	IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
	100	IF(lwp) WRITE(numout,*) '~~~~~~~'
[1847]	101	!
[1975]	102	rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp) ! Brown & Campana parameter for semi-implicit hpg
[592]	103	ENDIF
[1870]	104	! ! set time step size (Euler/Leapfrog)
	105	IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt_t(:) = rdttra(:) ! at nit000 (Euler)
	106	ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt_t(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog)
	107	ENDIF
[592]	108
[1870]	109	! !== Update after tracer on domain lateral boundaries ==!
[1110]	110	!
[1870]	111	CALL lbc_lnk( ta, 'T', 1. ) ! local domain boundaries (T-point, unchanged sign)
[3]	112	CALL lbc_lnk( sa, 'T', 1. )
[1110]	113	!
[3]	114	#if defined key_obc
[1870]	115	CALL obc_tra( kt ) ! OBC open boundaries
[3]	116	#endif
[1110]	117	#if defined key_bdy
[1870]	118	CALL bdy_tra( kt ) ! BDY open boundaries
[1110]	119	#endif
[392]	120	#if defined key_agrif
[1870]	121	CALL Agrif_tra ! AGRIF zoom boundaries
[389]	122	#endif
[1438]	123
[1870]	124	IF( l_trdtra ) THEN ! trends computation: store now fields before applying the Asselin filter
[1110]	125	ztrdt(:,:,:) = tn(:,:,:)
	126	ztrds(:,:,:) = sn(:,:,:)
	127	ENDIF
	128
[1870]	129	! !== modifed Leap-Frog + Asselin filter (modified LF-RA) scheme ==!
[1438]	130	IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt ) ! variable volume level (vvl)
	131	ELSE ; CALL tra_nxt_fix( kt ) ! fixed volume level
	132	ENDIF
	133
	134	#if defined key_agrif
[1870]	135	IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! Update tracer at AGRIF zoom boundaries (children only)
[1438]	136	#endif
	137
[1870]	138	IF( l_trdtra ) THEN ! trends computation: trend of the Asselin filter (tb filtered - tb)/dt
[1110]	139	DO jk = 1, jpkm1
[1438]	140	zfact = 1.e0 / r2dt_t(jk)
	141	ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact
	142	ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact
[1110]	143	END DO
[1438]	144	CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt )
	145	END IF
	146
	147	! ! control print
	148	IF(ln_ctl) CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt - Tn: ', mask1=tmask, &
	149	& tab3d_2=sn, clinfo2= ' Sn: ', mask2=tmask )
	150	!
	151	END SUBROUTINE tra_nxt
	152
	153
	154	SUBROUTINE tra_nxt_fix( kt )
	155	!!----------------------------------------------------------------------
	156	!! * ROUTINE tra_nxt_fix *
	157	!!
	158	!! ** Purpose : fixed volume: apply the Asselin time filter and
	159	!! swap the tracer fields.
	160	!!
	161	!! ** Method : - Apply a Asselin time filter on now fields.
	162	!! - save in (ta,sa) an average over the three time levels
	163	!! which will be used to compute rdn and thus the semi-implicit
	164	!! hydrostatic pressure gradient (ln_dynhpg_imp = T)
	165	!! - swap tracer fields to prepare the next time_step.
[1847]	166	!! This can be summurized for temperature as:
[1975]	167	!! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T
[1847]	168	!! ztm = 0 otherwise
[1975]	169	!! with rbcp=1/4 * (1-atfp^4) / (1-atfp)
[1438]	170	!! tb = tn + atfp*[ tb - 2 tn + ta ]
	171	!! tn = ta
	172	!! ta = ztm (NB: reset to 0 after eos_bn2 call)
	173	!!
	174	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	175	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	176	!!----------------------------------------------------------------------
[1870]	177	INTEGER, INTENT(in) :: kt ! ocean time-step index
[1438]	178	!!
[1870]	179	INTEGER :: ji, jj, jk ! dummy loop indices
[2068]	180	REAL(wp) :: zt_d, zs_d ! temporary scalars
[1975]	181	REAL(wp) :: ztn, zsn ! - -
[1438]	182	!!----------------------------------------------------------------------
	183
	184	IF( kt == nit000 ) THEN
	185	IF(lwp) WRITE(numout,*)
	186	IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
	187	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	188	ENDIF
	189	!
[1975]	190	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step
	191	DO jk = 1, jpkm1 ! (only swap)
	192	tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta
	193	sn(:,:,jk) = sa(:,:,jk) ! sn <-- sa
	194	END DO
	195	ELSE ! General case (Leapfrog + Asselin filter)
	196	DO jk = 1, jpkm1
	197	DO jj = 1, jpj
	198	DO ji = 1, jpi
	199	IF( ln_dynhpg_imp ) THEN ! implicit hpg: keep tn, sn in memory
	200	ztn = tn(ji,jj,jk)
	201	zsn = sn(ji,jj,jk)
	202	ENDIF
	203	! ! time laplacian on tracers
	204	! ! used for both Asselin and Brown & Campana filters
[2068]	205	zt_d = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk)
	206	zs_d = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk)
[1975]	207	!
	208	! ! swap of arrays
[2068]	209	tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_d ! tb <-- tn filtered
	210	sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_d ! sb <-- sn filtered
[1975]	211	tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta
	212	sn(ji,jj,jk) = sa(ji,jj,jk) ! sn <-- sa
	213	! ! semi imlicit hpg computation (Brown & Campana)
	214	IF( ln_dynhpg_imp ) THEN
[2068]	215	ta(ji,jj,jk) = ztn + rbcp * zt_d ! ta <-- Brown & Campana average
	216	sa(ji,jj,jk) = zsn + rbcp * zs_d ! sa <-- Brown & Campana average
[1975]	217	ENDIF
[3]	218	END DO
[1110]	219	END DO
[1975]	220	END DO
[1110]	221	ENDIF
[1438]	222	!
	223	END SUBROUTINE tra_nxt_fix
[3]	224
[1110]	225
[1438]	226	SUBROUTINE tra_nxt_vvl( kt )
	227	!!----------------------------------------------------------------------
	228	!! * ROUTINE tra_nxt_vvl *
	229	!!
	230	!! ** Purpose : Time varying volume: apply the Asselin time filter
	231	!! and swap the tracer fields.
	232	!!
	233	!! ** Method : - Apply a thickness weighted Asselin time filter on now fields.
	234	!! - save in (ta,sa) a thickness weighted average over the three
	235	!! time levels which will be used to compute rdn and thus the semi-
	236	!! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
	237	!! - swap tracer fields to prepare the next time_step.
[1847]	238	!! This can be summurized for temperature as:
[1975]	239	!! ztm = ( e3t_ntn + rbcp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] ) ln_dynhpg_imp = T
	240	!! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] )
[1847]	241	!! ztm = 0 otherwise
[1438]	242	!! tb = ( e3t_ntn + atfp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] )
	243	!! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] )
	244	!! tn = ta
	245	!! ta = zt (NB: reset to 0 after eos_bn2 call)
	246	!!
	247	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	248	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	249	!!----------------------------------------------------------------------
[1975]	250	INTEGER, INTENT(in) :: kt ! ocean time-step index
[1438]	251	!!
[1975]	252	INTEGER :: ji, jj, jk ! dummy loop indices
	253	REAL :: ze3t_a, ze3t_n, ze3t_b ! temporary scalars
	254	REAL :: ztc_a, ztc_n, ztc_b ! - -
	255	REAL :: zsc_a, zsc_n, zsc_b ! - -
[2068]	256	REAL :: ztc_f, zsc_f, ztc_d, zsc_d ! - -
	257	REAL :: ze3t_f, ze3t_d ! - -
[1975]	258	REAL :: zfact1, zfact2 ! - -
[1438]	259	!!----------------------------------------------------------------------
[1975]	260
[1438]	261	IF( kt == nit000 ) THEN
	262	IF(lwp) WRITE(numout,*)
	263	IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
	264	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	265	ENDIF
	266
[1975]	267	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step
	268	! ! (only swap)
	269	DO jk = 1, jpkm1 ! tn <-- ta
	270	tn(:,:,jk) = ta(:,:,jk) ! sn <-- sa
	271	sn(:,:,jk) = sa(:,:,jk)
	272	END DO
	273	! ! General case (Leapfrog + Asselin filter)
	274	ELSE ! apply filter on thickness weighted tracer and swap
	275	DO jk = 1, jpkm1
[2068]	276	zfact1 = atfp * rdttra(jk)
[1975]	277	zfact2 = zfact1 / rau0
	278	DO jj = 1, jpj
	279	DO ji = 1, jpi
	280	! ! scale factors at Before, now and after
	281	ze3t_b = fse3t_b(ji,jj,jk)
	282	ze3t_n = fse3t_n(ji,jj,jk)
	283	ze3t_a = fse3t_a(ji,jj,jk)
[2068]	284	ze3t_d = fse3t_d(ji,jj,jk)
[1975]	285	! ! tracer content at Before, now and after
	286	ztc_b = tb(ji,jj,jk) * ze3t_b ; zsc_b = sb(ji,jj,jk) * ze3t_b
	287	ztc_n = tn(ji,jj,jk) * ze3t_n ; zsc_n = sn(ji,jj,jk) * ze3t_n
	288	ztc_a = ta(ji,jj,jk) * ze3t_a ; zsc_a = sa(ji,jj,jk) * ze3t_a
	289	!
[2005]	290	! ! Time laplacian on tracer contents
[1975]	291	! ! used for both Asselin and Brown & Campana filters
[2068]	292	ztc_d = ztc_a - 2. * ztc_n + ztc_b
	293	zsc_d = zsc_a - 2. * zsc_n + zsc_b
[2005]	294	! ! Asselin Filter on thicknesses and tracer contents
[2068]	295	ze3t_f = ze3t_n + atfp * ze3t_d
	296	ztc_f = ztc_n + atfp * ztc_d
	297	zsc_f = zsc_n + atfp * zsc_d
[2005]	298	! ! Filter correction
[1975]	299	IF( jk == 1 ) THEN
[2091]	300	ze3t_f = ze3t_f - zfact2 * ( emp_b (ji,jj) - emp (ji,jj) )
	301	ztc_f = ztc_f - zfact1 * ( sbc_hc_n(ji,jj) - sbc_hc_b(ji,jj) )
	302	zsc_f = zsc_f - zfact1 * ( sbc_sc_n(ji,jj) - sbc_sc_b(ji,jj) )
[1975]	303	ENDIF
	304	IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN
[2091]	305	ztc_f = ztc_f - zfact1 * ( qsr_hc_n(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )
[1975]	306	ENDIF
[2068]	307	! swap of arrays
[1975]	308	ze3t_f = 1.e0 / ze3t_f
	309	tb(ji,jj,jk) = ztc_f * ze3t_f ! tb <-- tn filtered
	310	sb(ji,jj,jk) = zsc_f * ze3t_f ! sb <-- sn filtered
	311	tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta
	312	sn(ji,jj,jk) = sa(ji,jj,jk) ! sn <-- sa
	313	! ! semi imlicit hpg computation (Brown & Campana)
	314	IF( ln_dynhpg_imp ) THEN
[2068]	315	ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d )
	316	ta(ji,jj,jk) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average
	317	sa(ji,jj,jk) = ze3t_d * ( zsc_n + rbcp * zsc_d ) ! sa <-- Brown & Campana average
[1975]	318	ENDIF
[1438]	319	END DO
	320	END DO
[1975]	321	END DO
[1438]	322	ENDIF
[503]	323	!
[1438]	324	END SUBROUTINE tra_nxt_vvl
[3]	325
	326	!!======================================================================
	327	END MODULE tranxt

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

Download in other formats: