Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

tranxt.F90

Last change on this file was 6662, checked in by jgraham, 8 years ago
AMM15 is not running with CICE, so need to prevent unrealistically low temperatures.
File size: 18.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
[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
[4990]	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
[4990]	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)
[4990]	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
[4990]	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),
[4990]	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	!!
[6662]	95	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2528]	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
[6662]	110	! JGraham: Update after tracer to remove any temp < -2C (no sea ice in AMM15).
	111	DO jk = 1, jpk
	112	DO jj = 1, jpj
	113	DO ji = 1, jpi
	114	IF(tsa(ji,jj,jk,jp_tem) < -2.e0) THEN
	115	tsa(ji,jj,jk,jp_tem) = -2.e0
	116	ENDIF
	117	END DO
	118	END DO
	119	END DO
	120
[1110]	121	! Update after tracer on domain lateral boundaries
	122	!
[4230]	123	CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp ) ! local domain boundaries (T-point, unchanged sign)
	124	CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
[1110]	125	!
[2528]	126	#if defined key_bdy
[3294]	127	IF( lk_bdy ) CALL bdy_tra( kt ) ! BDY open boundaries
[1110]	128	#endif
[392]	129	#if defined key_agrif
[2528]	130	CALL Agrif_tra ! AGRIF zoom boundaries
[389]	131	#endif
[1438]	132
	133	! set time step size (Euler/Leapfrog)
[2715]	134	IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dtra(:) = rdttra(:) ! at nit000 (Euler)
[4230]	135	ELSEIF( kt <= nit000 + 1 ) THEN ; r2dtra(:) = 2._wp* rdttra(:) ! at nit000 or nit000+1 (Leapfrog)
[1438]	136	ENDIF
[3]	137
[1110]	138	! trends computation initialisation
[2528]	139	IF( l_trdtra ) THEN ! store now fields before applying the Asselin filter
[3294]	140	CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
	141	ztrdt(:,:,:) = tsn(:,:,:,jp_tem)
	142	ztrds(:,:,:) = tsn(:,:,:,jp_sal)
[4990]	143	IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend
	144	CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
	145	CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
	146	ENDIF
[1110]	147	ENDIF
	148
[2528]	149	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step (only swap)
	150	DO jn = 1, jpts
	151	DO jk = 1, jpkm1
	152	tsn(:,:,jk,jn) = tsa(:,:,jk,jn)
	153	END DO
	154	END DO
	155	ELSE ! Leap-Frog + Asselin filter time stepping
	156	!
[5385]	157	IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa, &
	158	& sbc_tsc, sbc_tsc_b, jpts ) ! variable volume level (vvl)
	159	ELSE ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! fixed volume level
[2528]	160	ENDIF
	161	ENDIF
[2715]	162	!
[1438]	163	#if defined key_agrif
	164	! Update tracer at AGRIF zoom boundaries
	165	IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! children only
	166	#endif
[2715]	167	!
[1438]	168	! trends computation
[2528]	169	IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt
[1110]	170	DO jk = 1, jpkm1
[4990]	171	zfact = 1._wp / r2dtra(jk)
[2528]	172	ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
	173	ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
[1110]	174	END DO
[4990]	175	CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
	176	CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
[3294]	177	CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
[1438]	178	END IF
[2715]	179	!
[1438]	180	! ! control print
[2528]	181	IF(ln_ctl) CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt - Tn: ', mask1=tmask, &
	182	& tab3d_2=tsn(:,:,:,jp_sal), clinfo2= ' Sn: ', mask2=tmask )
[1438]	183	!
[4990]	184	IF( nn_timing == 1 ) CALL timing_stop('tra_nxt')
[3294]	185	!
[1438]	186	END SUBROUTINE tra_nxt
	187
	188
[3294]	189	SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
[1438]	190	!!----------------------------------------------------------------------
	191	!! * ROUTINE tra_nxt_fix *
	192	!!
	193	!! ** Purpose : fixed volume: apply the Asselin time filter and
	194	!! swap the tracer fields.
	195	!!
	196	!! ** Method : - Apply a Asselin time filter on now fields.
	197	!! - save in (ta,sa) an average over the three time levels
	198	!! which will be used to compute rdn and thus the semi-implicit
	199	!! hydrostatic pressure gradient (ln_dynhpg_imp = T)
	200	!! - swap tracer fields to prepare the next time_step.
	201	!! This can be summurized for tempearture as:
[2528]	202	!! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T
	203	!! ztm = 0 otherwise
	204	!! with rbcp=1/4 * (1-atfp^4) / (1-atfp)
[1438]	205	!! tb = tn + atfp*[ tb - 2 tn + ta ]
[2528]	206	!! tn = ta
[1438]	207	!! ta = ztm (NB: reset to 0 after eos_bn2 call)
	208	!!
	209	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	210	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	211	!!----------------------------------------------------------------------
[2715]	212	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	213	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2715]	214	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	215	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	216	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields
	217	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields
	218	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
	219	!
[2528]	220	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	221	LOGICAL :: ll_tra_hpg ! local logical
	222	REAL(wp) :: ztn, ztd ! local scalars
[1438]	223	!!----------------------------------------------------------------------
	224
[3294]	225	IF( kt == kit000 ) THEN
[1438]	226	IF(lwp) WRITE(numout,*)
[3294]	227	IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
[1438]	228	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	229	ENDIF
	230	!
[2528]	231	IF( cdtype == 'TRA' ) THEN ; ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg
	232	ELSE ; ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg
	233	ENDIF
	234	!
	235	DO jn = 1, kjpt
[1438]	236	!
[2528]	237	DO jk = 1, jpkm1
	238	DO jj = 1, jpj
	239	DO ji = 1, jpi
	240	ztn = ptn(ji,jj,jk,jn)
	241	ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers
	242	!
	243	ptb(ji,jj,jk,jn) = ztn + atfp * ztd ! ptb <-- filtered ptn
	244	ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta
	245	!
	246	IF( ll_tra_hpg ) pta(ji,jj,jk,jn) = ztn + rbcp * ztd ! pta <-- Brown & Campana average
[3]	247	END DO
[2528]	248	END DO
	249	END DO
[1110]	250	!
[2528]	251	END DO
[1438]	252	!
	253	END SUBROUTINE tra_nxt_fix
[3]	254
[1110]	255
[5385]	256	SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
[1438]	257	!!----------------------------------------------------------------------
	258	!! * ROUTINE tra_nxt_vvl *
	259	!!
	260	!! ** Purpose : Time varying volume: apply the Asselin time filter
	261	!! and swap the tracer fields.
	262	!!
	263	!! ** Method : - Apply a thickness weighted Asselin time filter on now fields.
	264	!! - save in (ta,sa) a thickness weighted average over the three
	265	!! time levels which will be used to compute rdn and thus the semi-
	266	!! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
	267	!! - swap tracer fields to prepare the next time_step.
	268	!! This can be summurized for tempearture as:
[2528]	269	!! ztm = ( e3t_ntn + rbcp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] ) ln_dynhpg_imp = T
	270	!! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] )
	271	!! ztm = 0 otherwise
[1438]	272	!! tb = ( e3t_ntn + atfp[ e3t_btb - 2 e3t_ntn + e3t_a*ta ] )
	273	!! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] )
	274	!! tn = ta
	275	!! ta = zt (NB: reset to 0 after eos_bn2 call)
	276	!!
	277	!! ** Action : - (tb,sb) and (tn,sn) ready for the next time step
	278	!! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)
	279	!!----------------------------------------------------------------------
[5385]	280	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	281	INTEGER , INTENT(in ) :: kit000 ! first time step index
	282	REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! time-step
	283	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	284	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	285	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields
	286	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields
	287	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
	288	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc ! surface tracer content
	289	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc_b ! before surface tracer content
	290
[1438]	291	!!
[5467]	292	LOGICAL :: ll_tra_hpg, ll_traqsr, ll_rnf ! local logical
[2528]	293	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	294	REAL(wp) :: zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar
	295	REAL(wp) :: zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d ! - -
[1438]	296	!!----------------------------------------------------------------------
[3294]	297	!
	298	IF( kt == kit000 ) THEN
[1438]	299	IF(lwp) WRITE(numout,*)
[3294]	300	IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
[1438]	301	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	302	ENDIF
[2528]	303	!
	304	IF( cdtype == 'TRA' ) THEN
	305	ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg
	306	ll_traqsr = ln_traqsr ! active tracers case and solar penetration
[5467]	307	ll_rnf = ln_rnf ! active tracers case and river runoffs
[2528]	308	ELSE
	309	ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg
	310	ll_traqsr = .FALSE. ! active tracers case and NO solar penetration
[5467]	311	ll_rnf = .FALSE. ! passive tracers or NO river runoffs
[2528]	312	ENDIF
	313	!
	314	DO jn = 1, kjpt
	315	DO jk = 1, jpkm1
[5385]	316	zfact1 = atfp * p2dt(jk)
[2528]	317	zfact2 = zfact1 / rau0
	318	DO jj = 1, jpj
	319	DO ji = 1, jpi
	320	ze3t_b = fse3t_b(ji,jj,jk)
	321	ze3t_n = fse3t_n(ji,jj,jk)
	322	ze3t_a = fse3t_a(ji,jj,jk)
	323	! ! tracer content at Before, now and after
	324	ztc_b = ptb(ji,jj,jk,jn) * ze3t_b
	325	ztc_n = ptn(ji,jj,jk,jn) * ze3t_n
	326	ztc_a = pta(ji,jj,jk,jn) * ze3t_a
	327	!
	328	ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
	329	ztc_d = ztc_a - 2. * ztc_n + ztc_b
	330	!
	331	ze3t_f = ze3t_n + atfp * ze3t_d
	332	ztc_f = ztc_n + atfp * ztc_d
	333	!
[5385]	334	IF( jk == 1 ) THEN ! first level
[5467]	335	ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) )
[5385]	336	ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
[2528]	337	ENDIF
[5467]	338
[2528]	339	IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & ! solar penetration (temperature only)
	340	& ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )
[1438]	341
[5467]	342	IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) ) & ! river runoffs
	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
	346	ze3t_f = 1.e0 / ze3t_f
	347	ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered
	348	ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta
	349	!
	350	IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only)
	351	ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d )
	352	pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average
	353	ENDIF
[1438]	354	END DO
	355	END DO
[2528]	356	END DO
	357	!
	358	END DO
[503]	359	!
[1438]	360	END SUBROUTINE tra_nxt_vvl
[3]	361
	362	!!======================================================================
	363	END MODULE tranxt

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_amm15_test/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

Download in other formats: