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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 @ 2590

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

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

Download in other formats: