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.

trcsub.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90 @ 10963

Last change on this file since 10963 was 10963, checked in by acc, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert TOP routines in top-level TOP directory and all knock on effects of these conversions. SETTE tested.
Property svn:keywords set to `Id`
File size: 29.9 KB

Rev	Line
[2892]	1	MODULE trcsub
	2	!!======================================================================
	3	!! * MODULE trcsubstp *
[9019]	4	!! TOP : Averages physics variables for TOP substepping.
[2892]	5	!!======================================================================
	6	!! History : 1.0 ! 2011-10 (K. Edwards) Original
	7	!!----------------------------------------------------------------------
	8	#if defined key_top
	9	!!----------------------------------------------------------------------
[9019]	10	!! trc_sub : passive tracer system sub-stepping
[2892]	11	!!----------------------------------------------------------------------
[9019]	12	USE oce_trc ! ocean dynamics and active tracers variables
[2892]	13	USE trc
[9019]	14	USE trabbl ! bottom boundary layer
	15	USE zdf_oce
	16	USE domvvl
	17	USE divhor ! horizontal divergence
	18	USE sbcrnf , ONLY: h_rnf, nk_rnf ! River runoff
	19	USE bdy_oce , ONLY: ln_bdy, bdytmask ! BDY
	20	!
	21	USE prtctl_trc ! Print control for debbuging
	22	USE in_out_manager !
[2892]	23	USE iom
	24	USE lbclnk
	25	#if defined key_agrif
[9570]	26	USE agrif_oce_update
	27	USE agrif_oce_interp
[2892]	28	#endif
	29
	30	IMPLICIT NONE
	31
[9019]	32	PUBLIC trc_sub_stp ! called by trc_stp
	33	PUBLIC trc_sub_ini ! called by trc_ini to initialize substepping arrays.
	34	PUBLIC trc_sub_reset ! called by trc_stp to reset physics variables
	35	PUBLIC trc_sub_ssh ! called by trc_stp to reset physics variables
[2892]	36
[9019]	37	REAL(wp) :: r1_ndttrc ! = 1 / nn_dttrc
	38	REAL(wp) :: r1_ndttrcp1 ! = 1 / (nn_dttrc+1)
[2910]	39
[5836]	40
[9019]	41	!! averaged and temporary saved variables (needed when a larger passive tracer time-step is used)
	42	!! ----------------------------------------------------------------
	43	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_tm , un_temp !: i-horizontal velocity average [m/s]
	44	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vn_tm , vn_temp !: j-horizontal velocity average [m/s]
[10963]	45	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wn_temp !: hold current values of avt, uu(Kmm), vv(Kmm), ww
[9019]	46	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: tsn_tm , tsn_temp !: t/s average [m/s]
	47	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avs_tm , avs_temp !: vertical diffusivity coeff. at w-point [m2/s]
	48	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rhop_tm , rhop_temp !:
	49	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshn_tm , sshn_temp !: average ssh for the now step [m]
	50
	51	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf_tm , rnf_temp !: river runoff
	52	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_rnf_tm , h_rnf_temp !: depth in metres to the bottom of the relevant grid box
	53	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld_tm , hmld_temp !: mixed layer depth average [m]
	54	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr_i_tm , fr_i_temp !: average ice fraction [m/s]
	55	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_tm , emp_temp !: freshwater budget: volume flux [Kg/m2/s]
	56	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fmmflx_tm , fmmflx_temp !: freshwater budget: freezing/melting [Kg/m2/s]
	57	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_b_hold, emp_b_temp !: hold emp from the beginning of each sub-stepping[m]
	58	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr_tm , qsr_temp !: solar radiation average [m]
	59	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wndm_tm , wndm_temp !: 10m wind average [m]
	60	!
[10963]	61	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshb_hold !: hold ssh(Kbb) from the beginning of each sub-stepping[m]
[9019]	62	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshb_temp, ssha_temp
	63	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn_temp, rotn_temp
	64	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivb_temp, rotb_temp
	65	!
	66	! !!- bottom boundary layer param (ln_trabbl=T)
	67	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahu_bbl_tm, ahu_bbl_temp ! BBL diffusive i-coef.
	68	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahv_bbl_tm, ahv_bbl_temp ! BBL diffusive j-coef.
	69	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: utr_bbl_tm, utr_bbl_temp ! BBL u-advection
	70	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtr_bbl_tm, vtr_bbl_temp ! BBL v-advection
	71
	72	! !!- iso-neutral slopes (if l_ldfslp=T)
	73	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslp_temp, vslp_temp, wslpi_temp, wslpj_temp !: hold current values
	74	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslp_tm , vslp_tm , wslpi_tm , wslpj_tm !: time mean
	75
	76
[2892]	77	!!----------------------------------------------------------------------
[9598]	78	!! NEMO/TOP 4.0 , NEMO Consortium (2018)
[5215]	79	!! $Id$
[10068]	80	!! Software governed by the CeCILL license (see ./LICENSE)
[2892]	81	!!----------------------------------------------------------------------
	82	CONTAINS
	83
[10963]	84	SUBROUTINE trc_sub_stp( kt, Kbb, Kmm, Krhs )
[2892]	85	!!-------------------------------------------------------------------
	86	!! * ROUTINE trc_stp *
	87	!!
	88	!! ** Purpose : Average variables needed for sub-stepping passive tracers
	89	!!
	90	!! ** Method : Called every timestep to increment _tm (time mean) variables
	91	!! on TOP steps, calculate averages.
	92	!!-------------------------------------------------------------------
[10963]	93	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	94	INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time-level index
[9019]	95	!
	96	INTEGER :: ji, jj, jk ! dummy loop indices
	97	REAL(wp):: z1_ne3t, z1_ne3u, z1_ne3v, z1_ne3w ! local scalars
[2892]	98	!!-------------------------------------------------------------------
[3160]	99	!
[9124]	100	IF( ln_timing ) CALL timing_start('trc_sub_stp')
[3160]	101	!
	102	IF( kt == nit000 ) THEN
[2892]	103	IF(lwp) WRITE(numout,*)
	104	IF(lwp) WRITE(numout,*) 'trc_sub_stp : substepping of the passive tracers'
	105	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[2910]	106	!
[10963]	107	sshb_hold (:,:) = ssh (:,:,Kmm)
[2910]	108	emp_b_hold (:,:) = emp_b (:,:)
	109	!
	110	r1_ndttrc = 1._wp / REAL( nn_dttrc , wp )
	111	r1_ndttrcp1 = 1._wp / REAL( nn_dttrc + 1, wp )
[3160]	112	ENDIF
[2892]	113
[9019]	114	IF( MOD( kt , nn_dttrc ) /= 0 ) THEN
	115	!
[10963]	116	un_tm (:,:,:) = un_tm (:,:,:) + uu (:,:,:,Kmm) * e3u(:,:,:,Kmm)
	117	vn_tm (:,:,:) = vn_tm (:,:,:) + vv (:,:,:,Kmm) * e3v(:,:,:,Kmm)
	118	tsn_tm (:,:,:,jp_tem) = tsn_tm (:,:,:,jp_tem) + ts (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm)
	119	tsn_tm (:,:,:,jp_sal) = tsn_tm (:,:,:,jp_sal) + ts (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm)
	120	rhop_tm (:,:,:) = rhop_tm (:,:,:) + rhop (:,:,:) * e3t(:,:,:,Kmm)
	121	avs_tm (:,:,:) = avs_tm (:,:,:) + avs (:,:,:) * e3w(:,:,:,Kmm)
[5836]	122	IF( l_ldfslp ) THEN
	123	uslp_tm (:,:,:) = uslp_tm (:,:,:) + uslp (:,:,:)
	124	vslp_tm (:,:,:) = vslp_tm (:,:,:) + vslp (:,:,:)
	125	wslpi_tm(:,:,:) = wslpi_tm(:,:,:) + wslpi(:,:,:)
	126	wslpj_tm(:,:,:) = wslpj_tm(:,:,:) + wslpj(:,:,:)
	127	ENDIF
[9019]	128	IF( ln_trabbl ) THEN
	129	IF( nn_bbl_ldf == 1 ) THEN
	130	ahu_bbl_tm(:,:) = ahu_bbl_tm(:,:) + ahu_bbl(:,:)
	131	ahv_bbl_tm(:,:) = ahv_bbl_tm(:,:) + ahv_bbl(:,:)
	132	ENDIF
	133	IF( nn_bbl_adv == 1 ) THEN
	134	utr_bbl_tm(:,:) = utr_bbl_tm(:,:) + utr_bbl(:,:)
	135	vtr_bbl_tm(:,:) = vtr_bbl_tm(:,:) + vtr_bbl(:,:)
	136	ENDIF
	137	ENDIF
	138	!
[10963]	139	sshn_tm (:,:) = sshn_tm (:,:) + ssh (:,:,Kmm)
[9019]	140	rnf_tm (:,:) = rnf_tm (:,:) + rnf (:,:)
	141	h_rnf_tm (:,:) = h_rnf_tm (:,:) + h_rnf (:,:)
	142	hmld_tm (:,:) = hmld_tm (:,:) + hmld (:,:)
	143	fr_i_tm (:,:) = fr_i_tm (:,:) + fr_i (:,:)
	144	emp_tm (:,:) = emp_tm (:,:) + emp (:,:)
	145	fmmflx_tm(:,:) = fmmflx_tm(:,:) + fmmflx(:,:)
	146	qsr_tm (:,:) = qsr_tm (:,:) + qsr (:,:)
	147	wndm_tm (:,:) = wndm_tm (:,:) + wndm (:,:)
	148	!
[2910]	149	ELSE ! It is time to substep
[9019]	150	! 1. set temporary arrays to hold physics/dynamical variables
[10963]	151	un_temp (:,:,:) = uu (:,:,:,Kmm)
	152	vn_temp (:,:,:) = vv (:,:,:,Kmm)
	153	wn_temp (:,:,:) = ww (:,:,:)
	154	tsn_temp (:,:,:,:) = ts (:,:,:,:,Kmm)
[3192]	155	rhop_temp (:,:,:) = rhop (:,:,:)
	156	avs_temp (:,:,:) = avs (:,:,:)
[5836]	157	IF( l_ldfslp ) THEN
	158	uslp_temp (:,:,:) = uslp (:,:,:) ; wslpi_temp (:,:,:) = wslpi (:,:,:)
	159	vslp_temp (:,:,:) = vslp (:,:,:) ; wslpj_temp (:,:,:) = wslpj (:,:,:)
	160	ENDIF
[9019]	161	IF( ln_trabbl ) THEN
	162	IF( nn_bbl_ldf == 1 ) THEN
	163	ahu_bbl_temp(:,:) = ahu_bbl(:,:)
	164	ahv_bbl_temp(:,:) = ahv_bbl(:,:)
	165	ENDIF
	166	IF( nn_bbl_adv == 1 ) THEN
	167	utr_bbl_temp(:,:) = utr_bbl(:,:)
	168	vtr_bbl_temp(:,:) = vtr_bbl(:,:)
	169	ENDIF
	170	ENDIF
[10963]	171	sshn_temp (:,:) = ssh (:,:,Kmm)
	172	sshb_temp (:,:) = ssh (:,:,Kbb)
	173	ssha_temp (:,:) = ssh (:,:,Krhs)
[2910]	174	rnf_temp (:,:) = rnf (:,:)
	175	h_rnf_temp (:,:) = h_rnf (:,:)
	176	hmld_temp (:,:) = hmld (:,:)
[2944]	177	fr_i_temp (:,:) = fr_i (:,:)
[2910]	178	emp_temp (:,:) = emp (:,:)
	179	emp_b_temp (:,:) = emp_b (:,:)
[4148]	180	fmmflx_temp(:,:) = fmmflx(:,:)
[2944]	181	qsr_temp (:,:) = qsr (:,:)
	182	wndm_temp (:,:) = wndm (:,:)
[2971]	183	! ! Variables reset in trc_sub_ssh
[10963]	184	hdivn_temp (:,:,:) = hdiv (:,:,:)
[2910]	185	!
	186	! 2. Create averages and reassign variables
[10963]	187	un_tm (:,:,:) = un_tm (:,:,:) + uu (:,:,:,Kmm) * e3u(:,:,:,Kmm)
	188	vn_tm (:,:,:) = vn_tm (:,:,:) + vv (:,:,:,Kmm) * e3v(:,:,:,Kmm)
	189	tsn_tm (:,:,:,jp_tem) = tsn_tm (:,:,:,jp_tem) + ts (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm)
	190	tsn_tm (:,:,:,jp_sal) = tsn_tm (:,:,:,jp_sal) + ts (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm)
	191	rhop_tm (:,:,:) = rhop_tm (:,:,:) + rhop (:,:,:) * e3t(:,:,:,Kmm)
	192	avs_tm (:,:,:) = avs_tm (:,:,:) + avs (:,:,:) * e3w(:,:,:,Kmm)
[5836]	193	IF( l_ldfslp ) THEN
	194	uslp_tm (:,:,:) = uslp_tm (:,:,:) + uslp (:,:,:)
	195	vslp_tm (:,:,:) = vslp_tm (:,:,:) + vslp (:,:,:)
	196	wslpi_tm (:,:,:) = wslpi_tm(:,:,:) + wslpi(:,:,:)
	197	wslpj_tm (:,:,:) = wslpj_tm(:,:,:) + wslpj(:,:,:)
	198	ENDIF
[9019]	199	IF( ln_trabbl ) THEN
	200	IF( nn_bbl_ldf == 1 ) THEN
	201	ahu_bbl_tm(:,:) = ahu_bbl_tm(:,:) + ahu_bbl(:,:)
	202	ahv_bbl_tm(:,:) = ahv_bbl_tm(:,:) + ahv_bbl(:,:)
	203	ENDIF
	204	IF( nn_bbl_adv == 1 ) THEN
	205	utr_bbl_tm(:,:) = utr_bbl_tm(:,:) + utr_bbl(:,:)
	206	vtr_bbl_tm(:,:) = vtr_bbl_tm(:,:) + vtr_bbl(:,:)
	207	ENDIF
	208	ENDIF
[10963]	209	sshn_tm (:,:) = sshn_tm (:,:) + ssh (:,:,Kmm)
[2910]	210	rnf_tm (:,:) = rnf_tm (:,:) + rnf (:,:)
	211	h_rnf_tm (:,:) = h_rnf_tm (:,:) + h_rnf (:,:)
	212	hmld_tm (:,:) = hmld_tm (:,:) + hmld (:,:)
	213	fr_i_tm (:,:) = fr_i_tm (:,:) + fr_i (:,:)
	214	emp_tm (:,:) = emp_tm (:,:) + emp (:,:)
[4148]	215	fmmflx_tm(:,:) = fmmflx_tm (:,:) + fmmflx(:,:)
[2910]	216	qsr_tm (:,:) = qsr_tm (:,:) + qsr (:,:)
	217	wndm_tm (:,:) = wndm_tm (:,:) + wndm (:,:)
	218	!
[10963]	219	ssh (:,:,Kmm) = sshn_tm (:,:) * r1_ndttrcp1
	220	ssh (:,:,Kbb) = sshb_hold (:,:)
[2910]	221	rnf (:,:) = rnf_tm (:,:) * r1_ndttrcp1
	222	h_rnf (:,:) = h_rnf_tm (:,:) * r1_ndttrcp1
	223	hmld (:,:) = hmld_tm (:,:) * r1_ndttrcp1
[4611]	224	! variables that are initialized after averages
[10963]	225	emp_b (:,:) = emp_b_hold (:,:)
[2910]	226	IF( kt == nittrc000 ) THEN
	227	wndm (:,:) = wndm_tm (:,:) * r1_ndttrc
	228	qsr (:,:) = qsr_tm (:,:) * r1_ndttrc
	229	emp (:,:) = emp_tm (:,:) * r1_ndttrc
[4148]	230	fmmflx(:,:) = fmmflx_tm (:,:) * r1_ndttrc
[2944]	231	fr_i (:,:) = fr_i_tm (:,:) * r1_ndttrc
[9019]	232	IF( ln_trabbl ) THEN
	233	IF( nn_bbl_ldf == 1 ) THEN
	234	ahu_bbl(:,:) = ahu_bbl_tm (:,:) * r1_ndttrc
	235	ahv_bbl(:,:) = ahv_bbl_tm (:,:) * r1_ndttrc
	236	ENDIF
	237	IF( nn_bbl_adv == 1 ) THEN
	238	utr_bbl(:,:) = utr_bbl_tm (:,:) * r1_ndttrc
	239	vtr_bbl(:,:) = vtr_bbl_tm (:,:) * r1_ndttrc
	240	ENDIF
[2971]	241	ENDIF
[2910]	242	ELSE
	243	wndm (:,:) = wndm_tm (:,:) * r1_ndttrcp1
	244	qsr (:,:) = qsr_tm (:,:) * r1_ndttrcp1
	245	emp (:,:) = emp_tm (:,:) * r1_ndttrcp1
[4148]	246	fmmflx(:,:) = fmmflx_tm (:,:) * r1_ndttrcp1
[2944]	247	fr_i (:,:) = fr_i_tm (:,:) * r1_ndttrcp1
[9019]	248	IF( ln_trabbl ) THEN
	249	IF( nn_bbl_ldf == 1 ) THEN
	250	ahu_bbl(:,:) = ahu_bbl_tm (:,:) * r1_ndttrcp1
	251	ahv_bbl(:,:) = ahv_bbl_tm (:,:) * r1_ndttrcp1
	252	ENDIF
	253	IF( nn_bbl_adv == 1 ) THEN
	254	utr_bbl(:,:) = utr_bbl_tm (:,:) * r1_ndttrcp1
	255	vtr_bbl(:,:) = vtr_bbl_tm (:,:) * r1_ndttrcp1
	256	ENDIF
[2971]	257	ENDIF
[2910]	258	ENDIF
	259	!
[2892]	260	DO jk = 1, jpk
	261	DO jj = 1, jpj
	262	DO ji = 1, jpi
[10963]	263	z1_ne3t = r1_ndttrcp1 / e3t(ji,jj,jk,Kmm)
	264	z1_ne3u = r1_ndttrcp1 / e3u(ji,jj,jk,Kmm)
	265	z1_ne3v = r1_ndttrcp1 / e3v(ji,jj,jk,Kmm)
	266	z1_ne3w = r1_ndttrcp1 / e3w(ji,jj,jk,Kmm)
[2910]	267	!
[10963]	268	uu (ji,jj,jk,Kmm) = un_tm (ji,jj,jk) * z1_ne3u
	269	vv (ji,jj,jk,Kmm) = vn_tm (ji,jj,jk) * z1_ne3v
	270	ts (ji,jj,jk,jp_tem,Kmm) = tsn_tm (ji,jj,jk,jp_tem) * z1_ne3t
	271	ts (ji,jj,jk,jp_sal,Kmm) = tsn_tm (ji,jj,jk,jp_sal) * z1_ne3t
	272	rhop (ji,jj,jk) = rhop_tm (ji,jj,jk) * z1_ne3t
[9019]	273	!!gm : BUG ==>> for avs I don't understand the division by e3w
[10963]	274	avs (ji,jj,jk) = avs_tm (ji,jj,jk) * z1_ne3w
[5836]	275	END DO
	276	END DO
	277	END DO
	278	IF( l_ldfslp ) THEN
	279	wslpi(:,:,:) = wslpi_tm(:,:,:)
	280	wslpj(:,:,:) = wslpj_tm(:,:,:)
	281	uslp (:,:,:) = uslp_tm (:,:,:)
	282	vslp (:,:,:) = vslp_tm (:,:,:)
	283	ENDIF
[2944]	284	!
[10963]	285	CALL trc_sub_ssh( kt, Kbb, Kmm, Krhs ) ! after ssh & vertical velocity
[2944]	286	!
[2892]	287	ENDIF
[3160]	288	!
[9124]	289	IF( ln_timing ) CALL timing_stop('trc_sub_stp')
[3160]	290	!
[2892]	291	END SUBROUTINE trc_sub_stp
	292
[5836]	293
[10963]	294	SUBROUTINE trc_sub_ini( Kmm )
[2892]	295	!!-------------------------------------------------------------------
	296	!! * ROUTINE trc_sub_ini *
	297	!!
	298	!! ** Purpose : Initialize variables needed for sub-stepping passive tracers
	299	!!
	300	!! ** Method :
	301	!! Compute the averages for sub-stepping
	302	!!-------------------------------------------------------------------
	303	INTEGER :: ierr
[10963]	304	INTEGER, INTENT( in ) :: Kmm ! time level index
[2892]	305	!!-------------------------------------------------------------------
[3160]	306	!
[2892]	307	IF(lwp) WRITE(numout,*)
	308	IF(lwp) WRITE(numout,*) 'trc_sub_ini : initial set up of the passive tracers substepping'
	309	IF(lwp) WRITE(numout,*) '~~~~~~~'
	310
	311	ierr = trc_sub_alloc ()
[10425]	312	CALL mpp_sum( 'trcsub', ierr )
[2892]	313	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
	314
[10963]	315	un_tm (:,:,:) = uu (:,:,:,Kmm) * e3u(:,:,:,Kmm)
	316	vn_tm (:,:,:) = vv (:,:,:,Kmm) * e3v(:,:,:,Kmm)
	317	tsn_tm (:,:,:,jp_tem) = ts (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm)
	318	tsn_tm (:,:,:,jp_sal) = ts (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm)
	319	rhop_tm (:,:,:) = rhop (:,:,:) * e3t(:,:,:,Kmm)
[5836]	320	!!gm : BUG? ==>> for avt & avs I don't understand the division by e3w
[10963]	321	avs_tm (:,:,:) = avs (:,:,:) * e3w(:,:,:,Kmm)
[5836]	322	IF( l_ldfslp ) THEN
	323	wslpi_tm(:,:,:) = wslpi(:,:,:)
	324	wslpj_tm(:,:,:) = wslpj(:,:,:)
	325	uslp_tm (:,:,:) = uslp (:,:,:)
	326	vslp_tm (:,:,:) = vslp (:,:,:)
	327	ENDIF
[10963]	328	sshn_tm (:,:) = ssh (:,:,Kmm)
[2944]	329	rnf_tm (:,:) = rnf (:,:)
	330	h_rnf_tm (:,:) = h_rnf (:,:)
	331	hmld_tm (:,:) = hmld (:,:)
[2892]	332
[2910]	333	! Physics variables that are set after initialization:
[9019]	334	fr_i_tm (:,:) = 0._wp
	335	emp_tm (:,:) = 0._wp
[10963]	336	fmmflx_tm(:,:) = 0._wp
[9019]	337	qsr_tm (:,:) = 0._wp
	338	wndm_tm (:,:) = 0._wp
	339	IF( ln_trabbl ) THEN
	340	IF( nn_bbl_ldf == 1 ) THEN
	341	ahu_bbl_tm(:,:) = 0._wp
	342	ahv_bbl_tm(:,:) = 0._wp
	343	ENDIF
	344	IF( nn_bbl_adv == 1 ) THEN
	345	utr_bbl_tm(:,:) = 0._wp
	346	vtr_bbl_tm(:,:) = 0._wp
	347	ENDIF
[2971]	348	ENDIF
[2910]	349	!
[2892]	350	END SUBROUTINE trc_sub_ini
	351
[9124]	352
[10963]	353	SUBROUTINE trc_sub_reset( kt, Kbb, Kmm, Krhs )
[2892]	354	!!-------------------------------------------------------------------
	355	!! * ROUTINE trc_sub_reset *
	356	!!
	357	!! ** Purpose : Reset physics variables averaged for substepping
	358	!!
	359	!! ** Method :
	360	!! Compute the averages for sub-stepping
	361	!!-------------------------------------------------------------------
[10963]	362	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	363	INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs ! time level index
	364	INTEGER :: jk ! dummy loop indices
[2892]	365	!!-------------------------------------------------------------------
[3160]	366	!
[9124]	367	IF( ln_timing ) CALL timing_start('trc_sub_reset')
[3160]	368	!
[2910]	369	! restore physics variables
[10963]	370	uu (:,:,:,Kmm) = un_temp (:,:,:)
	371	vv (:,:,:,Kmm) = vn_temp (:,:,:)
	372	ww (:,:,:) = wn_temp (:,:,:)
	373	ts (:,:,:,:,Kmm) = tsn_temp (:,:,:,:)
	374	rhop (:,:,:) = rhop_temp(:,:,:)
	375	avs (:,:,:) = avs_temp (:,:,:)
[5836]	376	IF( l_ldfslp ) THEN
[10963]	377	wslpi (:,:,:) = wslpi_temp (:,:,:)
	378	wslpj (:,:,:) = wslpj_temp (:,:,:)
	379	uslp (:,:,:) = uslp_temp (:,:,:)
	380	vslp (:,:,:) = vslp_temp (:,:,:)
[5836]	381	ENDIF
[10963]	382	ssh (:,:,Kmm) = sshn_temp (:,:)
	383	ssh (:,:,Kbb) = sshb_temp (:,:)
	384	ssh (:,:,Krhs) = ssha_temp (:,:)
	385	rnf (:,:) = rnf_temp (:,:)
	386	h_rnf (:,:) = h_rnf_temp (:,:)
[2910]	387	!
	388	hmld (:,:) = hmld_temp (:,:)
	389	fr_i (:,:) = fr_i_temp (:,:)
	390	emp (:,:) = emp_temp (:,:)
[4148]	391	fmmflx(:,:) = fmmflx_temp(:,:)
[2910]	392	emp_b (:,:) = emp_b_temp (:,:)
[2944]	393	qsr (:,:) = qsr_temp (:,:)
	394	wndm (:,:) = wndm_temp (:,:)
[9019]	395	IF( ln_trabbl ) THEN
	396	IF( nn_bbl_ldf == 1 ) THEN
	397	ahu_bbl(:,:) = ahu_bbl_temp(:,:)
	398	ahv_bbl(:,:) = ahv_bbl_temp(:,:)
	399	ENDIF
	400	IF( nn_bbl_adv == 1 ) THEN
	401	utr_bbl(:,:) = utr_bbl_temp(:,:)
	402	vtr_bbl(:,:) = vtr_bbl_temp(:,:)
	403	ENDIF
[2971]	404	ENDIF
[2910]	405	!
[10963]	406	hdiv (:,:,:) = hdivn_temp (:,:,:)
[2910]	407	!
	408	! Start new averages
[10963]	409	un_tm (:,:,:) = uu (:,:,:,Kmm) * e3u(:,:,:,Kmm)
	410	vn_tm (:,:,:) = vv (:,:,:,Kmm) * e3v(:,:,:,Kmm)
	411	tsn_tm (:,:,:,jp_tem) = ts (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm)
	412	tsn_tm (:,:,:,jp_sal) = ts (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm)
	413	rhop_tm (:,:,:) = rhop (:,:,:) * e3t(:,:,:,Kmm)
	414	avs_tm (:,:,:) = avs (:,:,:) * e3w(:,:,:,Kmm)
[5836]	415	IF( l_ldfslp ) THEN
	416	uslp_tm (:,:,:) = uslp (:,:,:)
	417	vslp_tm (:,:,:) = vslp (:,:,:)
[4611]	418	wslpi_tm(:,:,:) = wslpi(:,:,:)
	419	wslpj_tm(:,:,:) = wslpj(:,:,:)
[5836]	420	ENDIF
[2910]	421	!
[10963]	422	sshb_hold (:,:) = ssh (:,:,Kmm )
[2944]	423	emp_b_hold (:,:) = emp (:,:)
[10963]	424	sshn_tm (:,:) = ssh (:,:,Kmm)
[2910]	425	rnf_tm (:,:) = rnf (:,:)
	426	h_rnf_tm (:,:) = h_rnf (:,:)
	427	hmld_tm (:,:) = hmld (:,:)
	428	fr_i_tm (:,:) = fr_i (:,:)
	429	emp_tm (:,:) = emp (:,:)
[4148]	430	fmmflx_tm (:,:) = fmmflx(:,:)
[2910]	431	qsr_tm (:,:) = qsr (:,:)
	432	wndm_tm (:,:) = wndm (:,:)
[9019]	433	IF( ln_trabbl ) THEN
	434	IF( nn_bbl_ldf == 1 ) THEN
	435	ahu_bbl_tm(:,:) = ahu_bbl(:,:)
	436	ahv_bbl_tm(:,:) = ahv_bbl(:,:)
	437	ENDIF
	438	IF( nn_bbl_adv == 1 ) THEN
	439	utr_bbl_tm(:,:) = utr_bbl(:,:)
	440	vtr_bbl_tm(:,:) = vtr_bbl(:,:)
	441	ENDIF
[2971]	442	ENDIF
[2910]	443	!
	444	!
[9124]	445	IF( ln_timing ) CALL timing_stop('trc_sub_reset')
[3160]	446	!
[2910]	447	END SUBROUTINE trc_sub_reset
[2892]	448
	449
[10963]	450	SUBROUTINE trc_sub_ssh( kt, Kbb, Kmm, Krhs )
[2892]	451	!!----------------------------------------------------------------------
	452	!! * ROUTINE trc_sub_ssh *
	453	!!
[10963]	454	!! ** Purpose : compute the after ssh (ssh(:,:,Krhs)), the now vertical velocity
[6140]	455	!! and update the now vertical coordinate (ln_linssh=F).
[2892]	456	!!
	457	!! ** Method : - Using the incompressibility hypothesis, the vertical
	458	!! velocity is computed by integrating the horizontal divergence
	459	!! from the bottom to the surface minus the scale factor evolution.
	460	!! The boundary conditions are w=0 at the bottom (no flux) and.
	461	!!
[10963]	462	!! ** action : ssh(:,:,Krhs) : after sea surface height
	463	!! ww : now vertical velocity
[6140]	464	!! sshu_a, sshv_a, sshf_a : after sea surface height (ln_linssh=F)
[2892]	465	!!
	466	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	467	!!----------------------------------------------------------------------
[10963]	468	INTEGER, INTENT(in) :: kt ! time step
	469	INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! ocean time-level index
[2892]	470	!
	471	INTEGER :: ji, jj, jk ! dummy loop indices
	472	REAL(wp) :: zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0 ! local scalars
[9125]	473	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
[3147]	474	!!---------------------------------------------------------------------
[3160]	475	!
[9124]	476	IF( ln_timing ) CALL timing_start('trc_sub_ssh')
[3160]	477	!
[2892]	478
	479	IF( kt == nittrc000 ) THEN
	480	!
	481	IF(lwp) WRITE(numout,*)
	482	IF(lwp) WRITE(numout,*) 'trc_sub_ssh : after sea surface height and now vertical velocity '
	483	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
	484	!
[10963]	485	ww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
[2892]	486	!
	487	ENDIF
	488	!
[10963]	489	!!gm BUG here ! hdiv will include the runoff divergence at the wrong timestep !!!!
[10922]	490	CALL div_hor( kt, Kmm ) ! Horizontal divergence & Relative vorticity
[2892]	491	!
	492	z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog)
	493	IF( neuler == 0 .AND. kt == nittrc000 ) z2dt = rdt
	494
	495	! !------------------------------!
	496	! ! After Sea Surface Height !
	497	! !------------------------------!
	498	zhdiv(:,:) = 0._wp
	499	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
[10963]	500	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
[2892]	501	END DO
	502	! ! Sea surface elevation time stepping
	503	! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
	504	! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
	505	z1_rau0 = 0.5 / rau0
[10963]	506	ssh(:,:,Krhs) = ( ssh(:,:,Kbb) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1)
[7646]	507
	508	IF( .NOT.ln_dynspg_ts ) THEN
[4611]	509	! These lines are not necessary with time splitting since
	510	! boundary condition on sea level is set during ts loop
[2892]	511	#if defined key_agrif
	512	CALL agrif_ssh( kt )
	513	#endif
[7646]	514	IF( ln_bdy ) THEN
[10963]	515	ssh(:,:,Krhs) = ssh(:,:,Krhs) * bdytmask(:,:)
	516	CALL lbc_lnk( 'trcsub', ssh(:,:,Krhs), 'T', 1. )
[7646]	517	ENDIF
	518	ENDIF
[6140]	519	!
[2892]	520	! !------------------------------!
	521	! ! Now Vertical Velocity !
	522	! !------------------------------!
	523	z1_2dt = 1.e0 / z2dt
	524	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
[6140]	525	! - ML - need 3 lines here because replacement of e3t by its expression yields too long lines otherwise
[10963]	526	ww(:,:,jk) = ww(:,:,jk+1) - e3t(:,:,jk,Kmm ) * hdiv(:,:,jk) &
	527	& - ( e3t(:,:,jk,Krhs) - e3t (:,:,jk,Kbb) ) &
[2892]	528	& * tmask(:,:,jk) * z1_2dt
[10963]	529	IF( ln_bdy ) ww(:,:,jk) = ww(:,:,jk) * bdytmask(:,:)
[2892]	530	END DO
	531	!
[9124]	532	IF( ln_timing ) CALL timing_stop('trc_sub_ssh')
[3160]	533	!
[2892]	534	END SUBROUTINE trc_sub_ssh
	535
[6140]	536
[2892]	537	INTEGER FUNCTION trc_sub_alloc()
	538	!!-------------------------------------------------------------------
	539	!! * ROUTINE trc_sub_alloc *
	540	!!-------------------------------------------------------------------
[10425]	541	USE lib_mpp, ONLY: ctl_stop
[9019]	542	INTEGER :: ierr(3)
[2892]	543	!!-------------------------------------------------------------------
	544	!
[9019]	545	ierr(:) = 0
[5836]	546	!
[9019]	547	ALLOCATE( un_temp(jpi,jpj,jpk) , vn_temp(jpi,jpj,jpk) , &
	548	& wn_temp(jpi,jpj,jpk) , &
	549	& rhop_temp(jpi,jpj,jpk) , rhop_tm(jpi,jpj,jpk) , &
	550	& sshn_temp(jpi,jpj) , sshb_temp(jpi,jpj) , &
	551	& ssha_temp(jpi,jpj) , &
	552	& rnf_temp(jpi,jpj) , h_rnf_temp(jpi,jpj) , &
	553	& tsn_temp(jpi,jpj,jpk,2) , emp_b_temp(jpi,jpj) , &
	554	& emp_temp(jpi,jpj) , fmmflx_temp(jpi,jpj) , &
	555	& hmld_temp(jpi,jpj) , qsr_temp(jpi,jpj) , &
	556	& fr_i_temp(jpi,jpj) , fr_i_tm(jpi,jpj) , &
	557	& wndm_temp(jpi,jpj) , wndm_tm(jpi,jpj) , &
	558	& avs_tm(jpi,jpj,jpk) , avs_temp(jpi,jpj,jpk) , &
	559	& hdivn_temp(jpi,jpj,jpk) , hdivb_temp(jpi,jpj,jpk), &
	560	& un_tm(jpi,jpj,jpk) , vn_tm(jpi,jpj,jpk) , &
	561	& sshn_tm(jpi,jpj) , sshb_hold(jpi,jpj) , &
	562	& tsn_tm(jpi,jpj,jpk,2) , &
	563	& emp_tm(jpi,jpj) , fmmflx_tm(jpi,jpj) , &
	564	& emp_b_hold(jpi,jpj) , &
	565	& hmld_tm(jpi,jpj) , qsr_tm(jpi,jpj) , &
	566	& rnf_tm(jpi,jpj) , h_rnf_tm(jpi,jpj) , STAT=ierr(1) )
[2892]	567	!
[5836]	568	IF( l_ldfslp ) THEN
[9019]	569	ALLOCATE( uslp_temp(jpi,jpj,jpk) , wslpi_temp(jpi,jpj,jpk), &
	570	& vslp_temp(jpi,jpj,jpk) , wslpj_temp(jpi,jpj,jpk), &
	571	& uslp_tm (jpi,jpj,jpk) , wslpi_tm (jpi,jpj,jpk), &
	572	& vslp_tm (jpi,jpj,jpk) , wslpj_tm (jpi,jpj,jpk), STAT=ierr(2) )
[5836]	573	ENDIF
[9019]	574	IF( ln_trabbl ) THEN
	575	ALLOCATE( ahu_bbl_temp(jpi,jpj) , utr_bbl_temp(jpi,jpj) , &
	576	& ahv_bbl_temp(jpi,jpj) , vtr_bbl_temp(jpi,jpj) , &
	577	& ahu_bbl_tm (jpi,jpj) , utr_bbl_tm (jpi,jpj) , &
	578	& ahv_bbl_tm (jpi,jpj) , vtr_bbl_tm (jpi,jpj) , STAT=ierr(3) )
	579	ENDIF
[5836]	580	!
[9019]	581	trc_sub_alloc = MAXVAL( ierr )
[5836]	582	!
[10425]	583	IF( trc_sub_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trc_sub_alloc: failed to allocate arrays' )
[9019]	584	!
[2892]	585	END FUNCTION trc_sub_alloc
	586
	587	#else
	588	!!----------------------------------------------------------------------
	589	!! Default key NO passive tracers
	590	!!----------------------------------------------------------------------
	591	CONTAINS
[10963]	592	SUBROUTINE trc_sub_stp( kt, Kbb, Kmm, Krhs ) ! Empty routine
	593	INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs ! time level indices
[2892]	594	WRITE(,) 'trc_sub_stp: You should not have seen this print! error?', kt
	595	END SUBROUTINE trc_sub_stp
[10963]	596	SUBROUTINE trc_sub_ini( Kmm ) ! Empty routine
	597	INTEGER, INTENT( in ) :: Kmm ! time level index
[2892]	598	WRITE(,) 'trc_sub_ini: You should not have seen this print! error?', kt
	599	END SUBROUTINE trc_sub_ini
	600	#endif
	601
	602	!!======================================================================
	603	END MODULE trcsub

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

Download in other formats: