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.

agrif_oce_update.F90 in NEMO/trunk/src/NST – NEMO

Context Navigation

source: NEMO/trunk/src/NST/agrif_oce_update.F90 @ 14083

Last change on this file since 14083 was 14064, checked in by ayoung, 4 years ago
Merging ticket #2506 into trunk.
Property svn:keywords set to `Id`
File size: 59.3 KB

Rev	Line
[12377]	1	#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES */
	2	#undef DECAL_FEEDBACK_2D /* SEPARATION of INTERFACES (Barotropic mode) */
	3	#undef VOL_REFLUX /* VOLUME REFLUXING*/
[5656]	4
[9570]	5	MODULE agrif_oce_update
[9019]	6	!!======================================================================
[9570]	7	!! * MODULE agrif_oce_interp *
[9019]	8	!! AGRIF: update package for the ocean dynamics (OPA)
	9	!!======================================================================
	10	!! History : 2.0 ! 2002-06 (L. Debreu) Original code
	11	!! 3.2 ! 2009-04 (R. Benshila)
	12	!! 3.6 ! 2014-09 (R. Benshila)
	13	!!----------------------------------------------------------------------
[7646]	14	#if defined key_agrif
[9019]	15	!!----------------------------------------------------------------------
	16	!! 'key_agrif' AGRIF zoom
	17	!!----------------------------------------------------------------------
[636]	18	USE par_oce
	19	USE oce
	20	USE dom_oce
[9019]	21	USE zdf_oce ! vertical physics: ocean variables
[782]	22	USE agrif_oce
[9019]	23	!
	24	USE in_out_manager ! I/O manager
	25	USE lib_mpp ! MPP library
[9031]	26	USE domvvl ! Need interpolation routines
[12377]	27	USE vremap ! Vertical remapping
[13216]	28	USE lbclnk
[14053]	29	#if defined key_qco
	30	USE domqco
	31	#endif
[636]	32	IMPLICIT NONE
	33	PRIVATE
[390]	34
[9055]	35	PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh
	36	PUBLIC Update_Scales
[9019]	37
[14053]	38	!! * Substitutions
	39	# include "domzgr_substitute.h90"
[1156]	40	!!----------------------------------------------------------------------
[9598]	41	!! NEMO/NST 4.0 , NEMO Consortium (2018)
[1156]	42	!! $Id$
[10068]	43	!! Software governed by the CeCILL license (see ./LICENSE)
[1156]	44	!!----------------------------------------------------------------------
[636]	45	CONTAINS
	46
[9031]	47	SUBROUTINE Agrif_Update_Tra( )
[9019]	48	!!----------------------------------------------------------------------
	49	!! * ROUTINE Agrif_Update_Tra *
	50	!!----------------------------------------------------------------------
[5656]	51	!
	52	IF (Agrif_Root()) RETURN
	53	!
[9031]	54	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update tracers from grid Number',Agrif_Fixed()
[636]	55
[12377]	56	#if defined key_vertical
	57	! Effect of this has to be carrefully checked
	58	! depending on what the nesting tools ensure for
	59	! volume conservation:
	60	Agrif_UseSpecialValueInUpdate = .FALSE.
	61	#else
[390]	62	Agrif_UseSpecialValueInUpdate = .TRUE.
[12377]	63	#endif
[9019]	64	Agrif_SpecialValueFineGrid = 0._wp
[5656]	65	!
	66	# if ! defined DECAL_FEEDBACK
[9031]	67	CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
	68	! near boundary update:
	69	! CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
[5656]	70	# else
[9031]	71	CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
	72	! near boundary update:
	73	! CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
[5656]	74	# endif
	75	!
[390]	76	Agrif_UseSpecialValueInUpdate = .FALSE.
[5656]	77	!
	78	!
[636]	79	END SUBROUTINE Agrif_Update_Tra
[390]	80
[9031]	81	SUBROUTINE Agrif_Update_Dyn( )
[9019]	82	!!----------------------------------------------------------------------
	83	!! * ROUTINE Agrif_Update_Dyn *
	84	!!----------------------------------------------------------------------
[5656]	85	!
	86	IF (Agrif_Root()) RETURN
	87	!
[9031]	88	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update momentum from grid Number',Agrif_Fixed()
[390]	89
[5656]	90	Agrif_UseSpecialValueInUpdate = .FALSE.
[13286]	91	Agrif_SpecialValueFineGrid = 0._wp
[13216]	92
	93	use_sign_north = .TRUE.
[13286]	94	sign_north = -1._wp
[13216]	95
[5656]	96	!
	97	# if ! defined DECAL_FEEDBACK
[9031]	98	CALL Agrif_Update_Variable(un_update_id,procname = updateU)
	99	CALL Agrif_Update_Variable(vn_update_id,procname = updateV)
	100	! near boundary update:
	101	! CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
	102	! CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)
[5656]	103	# else
[9031]	104	CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
	105	CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
	106	! near boundary update:
	107	! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
	108	! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
[5656]	109	# endif
[390]	110
[12377]	111	# if ! defined DECAL_FEEDBACK_2D
[5656]	112	CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
	113	CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)
	114	# else
	115	CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
	116	CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)
	117	# endif
[9031]	118	!
[12377]	119	# if ! defined DECAL_FEEDBACK_2D
[9031]	120	! Account for updated thicknesses at boundary edges
	121	IF (.NOT.ln_linssh) THEN
	122	! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)
	123	! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)
	124	ENDIF
	125	# endif
	126	!
[9019]	127	IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN
[4486]	128	! Update time integrated transports
[12377]	129	# if ! defined DECAL_FEEDBACK_2D
[9031]	130	CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
	131	CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
[5656]	132	# else
[9031]	133	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
	134	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
[5656]	135	# endif
	136	END IF
	137	!
[13216]	138	use_sign_north = .FALSE.
	139	!
[9031]	140	END SUBROUTINE Agrif_Update_Dyn
	141
	142	SUBROUTINE Agrif_Update_ssh( )
	143	!!---------------------------------------------
	144	!! * ROUTINE Agrif_Update_ssh *
	145	!!---------------------------------------------
	146	!
	147	IF (Agrif_Root()) RETURN
[5656]	148	!
	149	Agrif_UseSpecialValueInUpdate = .TRUE.
[13286]	150	Agrif_SpecialValueFineGrid = 0._wp
[12377]	151	# if ! defined DECAL_FEEDBACK_2D
[5656]	152	CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
	153	# else
	154	CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
	155	# endif
[9031]	156	!
[636]	157	Agrif_UseSpecialValueInUpdate = .FALSE.
[9031]	158	!
[9780]	159	# if defined VOL_REFLUX
[9031]	160	IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
[13216]	161	use_sign_north = .TRUE.
[13286]	162	sign_north = -1._wp
[9031]	163	! Refluxing on ssh:
[12377]	164	# if defined DECAL_FEEDBACK_2D
[9031]	165	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)
	166	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv)
[9780]	167	# else
	168	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu)
	169	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv)
	170	# endif
[13216]	171	use_sign_north = .FALSE.
[9031]	172	END IF
	173	# endif
	174	!
	175	END SUBROUTINE Agrif_Update_ssh
	176
	177
[9780]	178	SUBROUTINE Agrif_Update_Tke( )
[9031]	179	!!---------------------------------------------
	180	!! * ROUTINE Agrif_Update_Tke *
	181	!!---------------------------------------------
	182	!!
[5656]	183	!
[9031]	184	IF (Agrif_Root()) RETURN
	185	!
	186	Agrif_UseSpecialValueInUpdate = .TRUE.
	187	Agrif_SpecialValueFineGrid = 0.
	188
	189	CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN )
	190	CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
	191	CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
	192
	193	Agrif_UseSpecialValueInUpdate = .FALSE.
	194
	195	END SUBROUTINE Agrif_Update_Tke
	196
	197	SUBROUTINE Agrif_Update_vvl( )
	198	!!---------------------------------------------
	199	!! * ROUTINE Agrif_Update_vvl *
	200	!!---------------------------------------------
	201	!
	202	IF (Agrif_Root()) RETURN
	203	!
	204	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step()
	205	!
[14053]	206	#if ! defined key_qco
[9031]	207	Agrif_UseSpecialValueInUpdate = .TRUE.
	208	Agrif_SpecialValueFineGrid = 0.
	209	!
	210	! No interface separation here, update vertical grid at T points
	211	! everywhere over the overlapping regions (one account for refluxing in that case):
	212	CALL Agrif_Update_Variable(e3t_id, procname=updatee3t)
	213	!
	214	Agrif_UseSpecialValueInUpdate = .FALSE.
	215	!
	216	CALL Agrif_ChildGrid_To_ParentGrid()
	217	CALL dom_vvl_update_UVF
	218	CALL Agrif_ParentGrid_To_ChildGrid()
[14053]	219	#else
	220	CALL Agrif_ChildGrid_To_ParentGrid()
	221	CALL Agrif_Update_qco
	222	CALL Agrif_ParentGrid_To_ChildGrid()
	223	#endif
[9031]	224	!
	225	END SUBROUTINE Agrif_Update_vvl
	226
[14053]	227
	228	#if defined key_qco
	229	SUBROUTINE Agrif_Update_qco
	230	!!---------------------------------------------
	231	!! * ROUTINE dom_Update_qco *
	232	!!---------------------------------------------
	233	!
	234	! Save arrays prior update (needed for asselin correction)
	235	r3t(:,:,Krhs_a) = r3t(:,:,Kmm_a)
	236	r3u(:,:,Krhs_a) = r3u(:,:,Kmm_a)
	237	r3v(:,:,Krhs_a) = r3v(:,:,Kmm_a)
	238
	239	! Update r3x arrays from updated ssh
	240	CALL dom_qco_zgr( Kbb_a, Kmm_a )
	241	!
	242	END SUBROUTINE Agrif_Update_qco
	243	#endif
	244
	245
	246	#if ! defined key_qco
[9031]	247	SUBROUTINE dom_vvl_update_UVF
	248	!!---------------------------------------------
	249	!! * ROUTINE dom_vvl_update_UVF *
	250	!!---------------------------------------------
	251	!!
	252	INTEGER :: jk
	253	REAL(wp):: zcoef
	254	!!---------------------------------------------
	255	IF (lwp.AND.lk_agrif_debug) Write(,) 'Finalize e3 on grid Number', &
	256	& Agrif_Fixed(), 'Step', Agrif_Nb_Step()
	257
	258	! Save "old" scale factor (prior update) for subsequent asselin correction
	259	! of prognostic variables
	260	! -----------------------
	261	!
[12377]	262	e3u(:,:,:,Krhs_a) = e3u(:,:,:,Kmm_a)
	263	e3v(:,:,:,Krhs_a) = e3v(:,:,:,Kmm_a)
	264	! uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a)
	265	! vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a)
	266	hu(:,:,Krhs_a) = hu(:,:,Kmm_a)
	267	hv(:,:,Krhs_a) = hv(:,:,Kmm_a)
[9031]	268
	269	! 1) NOW fields
	270	!--------------
	271
	272	! Vertical scale factor interpolations
	273	! ------------------------------------
[12377]	274	CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3u(:,:,:,Kmm_a) , 'U' )
	275	CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3v(:,:,:,Kmm_a) , 'V' )
	276	CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3f(:,:,:) , 'F' )
[9031]	277
[12377]	278	CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3uw(:,:,:,Kmm_a), 'UW' )
	279	CALL dom_vvl_interpol( e3v(:,:,:,Kmm_a), e3vw(:,:,:,Kmm_a), 'VW' )
[9031]	280
	281	! Update total depths:
	282	! --------------------
[12377]	283	hu(:,:,Kmm_a) = 0._wp ! Ocean depth at U-points
	284	hv(:,:,Kmm_a) = 0._wp ! Ocean depth at V-points
[9031]	285	DO jk = 1, jpkm1
[12377]	286	hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk)
	287	hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk)
[9031]	288	END DO
	289	! ! Inverse of the local depth
[12377]	290	r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) )
	291	r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) )
[9031]	292
	293
	294	! 2) BEFORE fields:
	295	!------------------
[12489]	296	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN
[9031]	297	!
	298	! Vertical scale factor interpolations
	299	! ------------------------------------
[12377]	300	CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3u(:,:,:,Kbb_a), 'U' )
	301	CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3v(:,:,:,Kbb_a), 'V' )
[9031]	302
[12377]	303	CALL dom_vvl_interpol( e3u(:,:,:,Kbb_a), e3uw(:,:,:,Kbb_a), 'UW' )
	304	CALL dom_vvl_interpol( e3v(:,:,:,Kbb_a), e3vw(:,:,:,Kbb_a), 'VW' )
[9031]	305
	306	! Update total depths:
	307	! --------------------
[12377]	308	hu(:,:,Kbb_a) = 0._wp ! Ocean depth at U-points
	309	hv(:,:,Kbb_a) = 0._wp ! Ocean depth at V-points
[9031]	310	DO jk = 1, jpkm1
[12377]	311	hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk)
	312	hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk)
[9031]	313	END DO
	314	! ! Inverse of the local depth
[12377]	315	r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) )
	316	r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) )
[5656]	317	ENDIF
	318	!
[9031]	319	END SUBROUTINE dom_vvl_update_UVF
[14053]	320	#endif
[636]	321
[9031]	322	#if defined key_vertical
[6140]	323
[3294]	324	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
[9019]	325	!!----------------------------------------------------------------------
[636]	326	!! * ROUTINE updateT *
[9031]	327	!!---------------------------------------------
	328	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
	329	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	330	LOGICAL, INTENT(in) :: before
	331	!!
	332	INTEGER :: ji,jj,jk,jn
[12377]	333	INTEGER :: N_in, N_out
	334	REAL(wp) :: ztb, ztnu, ztno
[9031]	335	REAL(wp) :: h_in(k1:k2)
	336	REAL(wp) :: h_out(1:jpk)
[12377]	337	REAL(wp) :: tabin(k1:k2,1:jpts)
	338	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: tabres_child
[9031]	339	!!---------------------------------------------
[5656]	340	!
[9031]	341	IF (before) THEN
[12377]	342	!jc_alt
	343	! AGRIF_SpecialValue = -999._wp
[9031]	344	DO jn = n1,n2-1
	345	DO jk=k1,k2
	346	DO jj=j1,j2
	347	DO ji=i1,i2
[12377]	348	!jc_alt
	349	! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) &
	350	! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp
	351	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
[9031]	352	END DO
	353	END DO
	354	END DO
	355	END DO
	356	DO jk=k1,k2
	357	DO jj=j1,j2
	358	DO ji=i1,i2
[12377]	359	!jc_alt
	360	! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
	361	! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp
	362	tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
[9031]	363	END DO
	364	END DO
	365	END DO
	366	ELSE
[12377]	367	tabres_child(:,:,:,:) = 0._wp
[9031]	368	AGRIF_SpecialValue = 0._wp
	369	DO jj=j1,j2
	370	DO ji=i1,i2
	371	N_in = 0
	372	DO jk=k1,k2 !k2 = jpk of child grid
[12377]	373	! jc_alt
	374	! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT
	375	IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT
[9031]	376	N_in = N_in + 1
	377	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
	378	h_in(N_in) = tabres(ji,jj,jk,n2)
	379	ENDDO
	380	N_out = 0
	381	DO jk=1,jpk ! jpk of parent grid
[12377]	382	IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF
[9031]	383	N_out = N_out + 1
[12377]	384	h_out(N_out) = e3t(ji,jj,jk,Kmm_a)
[9031]	385	ENDDO
[12377]	386	IF (N_in*N_out > 0) THEN !Remove this?
	387	CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
[9031]	388	ENDIF
	389	ENDDO
	390	ENDDO
	391
[12489]	392	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
[9031]	393	! Add asselin part
[12377]	394	DO jn = 1,jpts
	395	DO jk = 1, jpkm1
	396	DO jj = j1, j2
	397	DO ji = i1, i2
	398	IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
	399	ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
	400	ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
	401	ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
[12489]	402	ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) &
[12377]	403	& * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
[9031]	404	ENDIF
[12377]	405	END DO
	406	END DO
	407	END DO
	408	END DO
[9031]	409	ENDIF
[12377]	410	DO jn = 1,jpts
	411	DO jk = 1, jpkm1
	412	DO jj = j1, j2
	413	DO ji = i1, i2
	414	IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
	415	ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn)
[9031]	416	END IF
	417	END DO
	418	END DO
	419	END DO
	420	END DO
[12377]	421	!
[12489]	422	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	423	ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a)
	424	ENDIF
[9031]	425	ENDIF
	426	!
	427	END SUBROUTINE updateTS
	428
	429	# else
	430
	431	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
	432	!!---------------------------------------------
	433	!! * ROUTINE updateT *
	434	!!---------------------------------------------
	435	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
	436	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	437	LOGICAL, INTENT(in) :: before
	438	!!
	439	INTEGER :: ji,jj,jk,jn
	440	REAL(wp) :: ztb, ztnu, ztno
	441	!!---------------------------------------------
[9019]	442	!
[9031]	443	IF (before) THEN
	444	DO jn = 1,jpts
	445	DO jk=k1,k2
	446	DO jj=j1,j2
	447	DO ji=i1,i2
	448	!> jc tmp
[12377]	449	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
	450	! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
[9031]	451	!< jc tmp
[3294]	452	END DO
[636]	453	END DO
	454	END DO
	455	END DO
	456	ELSE
[9031]	457	!> jc tmp
	458	DO jn = 1,jpts
	459	tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
	460	& * tmask(i1:i2,j1:j2,k1:k2)
	461	ENDDO
	462	!< jc tmp
[12489]	463	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
[5656]	464	! Add asselin part
[9031]	465	DO jn = 1,jpts
[9019]	466	DO jk = k1, k2
	467	DO jj = j1, j2
	468	DO ji = i1, i2
	469	IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
[12377]	470	ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
[9031]	471	ztnu = tabres(ji,jj,jk,jn)
[12377]	472	ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
[12489]	473	ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) &
[12377]	474	& * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
[4491]	475	ENDIF
[9019]	476	END DO
	477	END DO
	478	END DO
	479	END DO
[4491]	480	ENDIF
[9031]	481	DO jn = 1,jpts
[3294]	482	DO jk=k1,k2
	483	DO jj=j1,j2
	484	DO ji=i1,i2
[9019]	485	IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
[12377]	486	ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
[3294]	487	END IF
	488	END DO
[636]	489	END DO
	490	END DO
	491	END DO
[9031]	492	!
[12489]	493	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	494	ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a)
[9031]	495	ENDIF
	496	!
[636]	497	ENDIF
[5656]	498	!
[3294]	499	END SUBROUTINE updateTS
[390]	500
[9031]	501	# endif
[6140]	502
[9031]	503	# if defined key_vertical
	504
	505	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
[636]	506	!!---------------------------------------------
	507	!! * ROUTINE updateu *
	508	!!---------------------------------------------
[9031]	509	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	510	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	511	LOGICAL , INTENT(in ) :: before
[6140]	512	!
[9019]	513	INTEGER :: ji, jj, jk
[12377]	514	REAL(wp):: zrhoy, zub, zunu, zuno
[9031]	515	! VERTICAL REFINEMENT BEGIN
	516	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	517	REAL(wp) :: h_in(k1:k2)
	518	REAL(wp) :: h_out(1:jpk)
	519	INTEGER :: N_in, N_out
	520	REAL(wp) :: h_diff, excess, thick
	521	REAL(wp) :: tabin(k1:k2)
	522	! VERTICAL REFINEMENT END
[5656]	523	!!---------------------------------------------
	524	!
[6140]	525	IF( before ) THEN
[636]	526	zrhoy = Agrif_Rhoy()
[12377]	527	!jc_alt
	528	! AGRIF_SpecialValue = -999._wp
[9031]	529	DO jk=k1,k2
	530	DO jj=j1,j2
	531	DO ji=i1,i2
[12377]	532	!jc_alt
	533	! tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) &
	534	! & + (umask(ji,jj,jk)-1._wp)*999._wp
	535	tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)
	536	!jc_alt
	537	! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) &
	538	! & + (umask(ji,jj,jk)-1._wp)*999._wp
	539	tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
[9031]	540	END DO
	541	END DO
	542	END DO
	543	ELSE
	544	tabres_child(:,:,:) = 0.
	545	AGRIF_SpecialValue = 0._wp
	546	DO jj=j1,j2
	547	DO ji=i1,i2
	548	N_in = 0
	549	h_in(:) = 0._wp
	550	tabin(:) = 0._wp
	551	DO jk=k1,k2 !k2=jpk of child grid
[12377]	552	!jc_alt
	553	! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT
	554	IF( tabres(ji,jj,jk,2) == 0.) EXIT
[9031]	555	N_in = N_in + 1
	556	tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
	557	h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
	558	ENDDO
	559	N_out = 0
	560	DO jk=1,jpk
	561	IF (umask(ji,jj,jk) == 0) EXIT
	562	N_out = N_out + 1
[12377]	563	h_out(N_out) = e3u(ji,jj,jk,Kmm_a)
[9031]	564	ENDDO
	565	IF (N_in * N_out > 0) THEN
	566	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
[12377]	567	excess = 0._wp
[9031]	568	IF (h_diff < -1.e-4) THEN
	569	!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
	570	!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
	571	DO jk=N_in,1,-1
	572	thick = MIN(-1*h_diff, h_in(jk))
	573	excess = excess + tabin(jk)thicke2u(ji,jj)
	574	tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
	575	h_diff = h_diff + thick
	576	IF ( h_diff == 0) THEN
	577	N_in = jk
	578	h_in(jk) = h_in(jk) - thick
	579	EXIT
	580	ENDIF
	581	ENDDO
	582	ENDIF
[12377]	583	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
[9031]	584	tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
	585	ENDIF
	586	ENDDO
	587	ENDDO
[12377]	588	!
[9031]	589	DO jk=1,jpk
	590	DO jj=j1,j2
	591	DO ji=i1,i2
[12489]	592	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
[12377]	593	zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
	594	zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
	595	zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a)
[12489]	596	uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &
[12377]	597	& * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
[9031]	598	ENDIF
	599	!
[12377]	600	uu(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
[9031]	601	END DO
	602	END DO
	603	END DO
[12377]	604	!
[12489]	605	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	606	uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a)
	607	ENDIF
	608	!
[9031]	609	ENDIF
	610	!
	611	END SUBROUTINE updateu
	612
	613	#else
	614
	615	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
	616	!!---------------------------------------------
	617	!! * ROUTINE updateu *
	618	!!---------------------------------------------
	619	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	620	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	621	LOGICAL , INTENT(in ) :: before
	622	!
	623	INTEGER :: ji, jj, jk
	624	REAL(wp) :: zrhoy, zub, zunu, zuno
	625	!!---------------------------------------------
	626	!
	627	IF( before ) THEN
	628	zrhoy = Agrif_Rhoy()
[6140]	629	DO jk = k1, k2
[12377]	630	tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a)
[636]	631	END DO
	632	ELSE
[390]	633	DO jk=k1,k2
[636]	634	DO jj=j1,j2
	635	DO ji=i1,i2
[9031]	636	tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)
[4491]	637	!
[12489]	638	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
[12377]	639	zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
	640	zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
[9031]	641	zunu = tabres(ji,jj,jk,1)
[12489]	642	uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &
[12377]	643	& * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
[4491]	644	ENDIF
	645	!
[12377]	646	uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a)
[636]	647	END DO
	648	END DO
	649	END DO
[9031]	650	!
[12489]	651	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	652	uu(i1:i2,j1:j2,k1:k2,Kbb_a) = uu(i1:i2,j1:j2,k1:k2,Kmm_a)
[9031]	653	ENDIF
	654	!
[636]	655	ENDIF
[5656]	656	!
[636]	657	END SUBROUTINE updateu
[390]	658
[9031]	659	# endif
[6140]	660
[9031]	661	SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
	662	!!---------------------------------------------
	663	!! * ROUTINE correct_u_bdy *
	664	!!---------------------------------------------
	665	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	666	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	667	LOGICAL , INTENT(in ) :: before
	668	INTEGER , INTENT(in) :: nb, ndir
[5656]	669	!!
[9031]	670	LOGICAL :: western_side, eastern_side
	671	!
	672	INTEGER :: jj, jk
	673	REAL(wp) :: zcor
	674	!!---------------------------------------------
	675	!
	676	IF( .NOT.before ) THEN
	677	!
	678	western_side = (nb == 1).AND.(ndir == 1)
	679	eastern_side = (nb == 1).AND.(ndir == 2)
	680	!
	681	IF (western_side) THEN
	682	DO jj=j1,j2
[12377]	683	zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a)
	684	uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor
[9031]	685	DO jk=1,jpkm1
[12377]	686	uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk)
[9031]	687	END DO
	688	END DO
	689	ENDIF
	690	!
	691	IF (eastern_side) THEN
	692	DO jj=j1,j2
[12377]	693	zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a)
	694	uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor
[9031]	695	DO jk=1,jpkm1
[12377]	696	uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk)
[9031]	697	END DO
	698	END DO
	699	ENDIF
	700	!
	701	ENDIF
	702	!
	703	END SUBROUTINE correct_u_bdy
	704
	705	# if defined key_vertical
	706
	707	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
	708	!!---------------------------------------------
	709	!! * ROUTINE updatev *
	710	!!---------------------------------------------
	711	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	712	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	713	LOGICAL , INTENT(in ) :: before
	714	!
[9019]	715	INTEGER :: ji, jj, jk
[12377]	716	REAL(wp) :: zrhox, zvb, zvnu, zvno
[9031]	717	! VERTICAL REFINEMENT BEGIN
	718	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	719	REAL(wp) :: h_in(k1:k2)
	720	REAL(wp) :: h_out(1:jpk)
	721	INTEGER :: N_in, N_out
	722	REAL(wp) :: h_diff, excess, thick
	723	REAL(wp) :: tabin(k1:k2)
	724	! VERTICAL REFINEMENT END
	725	!!---------------------------------------------
[5656]	726	!
[9019]	727	IF( before ) THEN
[636]	728	zrhox = Agrif_Rhox()
[12377]	729	!jc_alt
	730	! AGRIF_SpecialValue = -999._wp
[9031]	731	DO jk=k1,k2
	732	DO jj=j1,j2
	733	DO ji=i1,i2
[12377]	734	!jc_alt
	735	! tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) &
	736	! & + (vmask(ji,jj,jk)-1._wp) * 999._wp
	737	tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a)
	738	!jc_alt
	739	! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) &
	740	! & + (vmask(ji,jj,jk)-1._wp) * 999._wp
	741	tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk)
[636]	742	END DO
	743	END DO
	744	END DO
	745	ELSE
[9031]	746	tabres_child(:,:,:) = 0.
	747	AGRIF_SpecialValue = 0._wp
	748	DO jj=j1,j2
	749	DO ji=i1,i2
	750	N_in = 0
	751	DO jk=k1,k2
[12377]	752	!jc_alt
	753	! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT
	754	IF (tabres(ji,jj,jk,2) == 0) EXIT
[9031]	755	N_in = N_in + 1
	756	tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
	757	h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
	758	ENDDO
	759	N_out = 0
	760	DO jk=1,jpk
	761	IF (vmask(ji,jj,jk) == 0) EXIT
	762	N_out = N_out + 1
[12377]	763	h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
[9031]	764	ENDDO
	765	IF (N_in * N_out > 0) THEN
	766	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
[12377]	767	excess = 0._wp
[9031]	768	IF (h_diff < -1.e-4) then
[12377]	769	!Even if bathy at T points match it's possible for the V points to be deeper in the child grid.
[9031]	770	!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
	771	DO jk=N_in,1,-1
	772	thick = MIN(-1*h_diff, h_in(jk))
	773	excess = excess + tabin(jk)thicke2u(ji,jj)
	774	tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
	775	h_diff = h_diff + thick
	776	IF ( h_diff == 0) THEN
	777	N_in = jk
	778	h_in(jk) = h_in(jk) - thick
	779	EXIT
	780	ENDIF
	781	ENDDO
	782	ENDIF
[12377]	783	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
[9031]	784	tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
	785	ENDIF
	786	ENDDO
	787	ENDDO
[12377]	788	!
	789	DO jk=1,jpkm1
[9031]	790	DO jj=j1,j2
	791	DO ji=i1,i2
[12489]	792	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
[12377]	793	zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
	794	zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
	795	zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a)
[12489]	796	vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &
[12377]	797	& * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
[4491]	798	ENDIF
	799	!
[12377]	800	vv(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
[636]	801	END DO
	802	END DO
	803	END DO
[12377]	804	!
[12489]	805	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	806	vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a)
	807	ENDIF
	808	!
[636]	809	ENDIF
[5656]	810	!
[636]	811	END SUBROUTINE updatev
[390]	812
[9031]	813	# else
[6140]	814
[9031]	815	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
	816	!!---------------------------------------------
	817	!! * ROUTINE updatev *
	818	!!---------------------------------------------
	819	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	820	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	821	LOGICAL , INTENT(in ) :: before
	822	!
	823	INTEGER :: ji, jj, jk
	824	REAL(wp) :: zrhox, zvb, zvnu, zvno
	825	!!---------------------------------------------
	826	!
	827	IF (before) THEN
	828	zrhox = Agrif_Rhox()
	829	DO jk=k1,k2
	830	DO jj=j1,j2
	831	DO ji=i1,i2
[12377]	832	tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
[9031]	833	END DO
	834	END DO
	835	END DO
	836	ELSE
	837	DO jk=k1,k2
	838	DO jj=j1,j2
	839	DO ji=i1,i2
	840	tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
	841	!
[12489]	842	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
[12377]	843	zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
	844	zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
[9031]	845	zvnu = tabres(ji,jj,jk,1)
[12489]	846	vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &
[12377]	847	& * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
[9031]	848	ENDIF
	849	!
[12377]	850	vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a)
[9031]	851	END DO
	852	END DO
	853	END DO
	854	!
[12489]	855	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	856	vv(i1:i2,j1:j2,k1:k2,Kbb_a) = vv(i1:i2,j1:j2,k1:k2,Kmm_a)
[9031]	857	ENDIF
	858	!
	859	ENDIF
	860	!
	861	END SUBROUTINE updatev
	862
	863	# endif
	864
	865	SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
	866	!!---------------------------------------------
[13216]	867	!! * ROUTINE correct_v_bdy *
[9031]	868	!!---------------------------------------------
	869	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	870	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	871	LOGICAL , INTENT(in ) :: before
	872	INTEGER , INTENT(in) :: nb, ndir
	873	!!
	874	LOGICAL :: southern_side, northern_side
	875	!
	876	INTEGER :: ji, jk
	877	REAL(wp) :: zcor
	878	!!---------------------------------------------
	879	!
	880	IF( .NOT.before ) THEN
	881	!
	882	southern_side = (nb == 2).AND.(ndir == 1)
	883	northern_side = (nb == 2).AND.(ndir == 2)
	884	!
	885	IF (southern_side) THEN
	886	DO ji=i1,i2
[12377]	887	zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a)
	888	vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor
[9031]	889	DO jk=1,jpkm1
[12377]	890	vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk)
[9031]	891	END DO
	892	END DO
	893	ENDIF
	894	!
	895	IF (northern_side) THEN
	896	DO ji=i1,i2
[12377]	897	zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a)
	898	vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor
[9031]	899	DO jk=1,jpkm1
[12377]	900	vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk)
[9031]	901	END DO
	902	END DO
	903	ENDIF
	904	!
	905	ENDIF
	906	!
	907	END SUBROUTINE correct_v_bdy
	908
	909
[636]	910	SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
[9019]	911	!!----------------------------------------------------------------------
	912	!! * ROUTINE updateu2d *
	913	!!----------------------------------------------------------------------
	914	INTEGER , INTENT(in ) :: i1, i2, j1, j2
	915	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	916	LOGICAL , INTENT(in ) :: before
[14064]	917	!!
	918	REAL(wp), DIMENSION(jpi,jpj) :: zpgu ! 2D workspace
[5656]	919	!!
[9031]	920	INTEGER :: ji, jj, jk
	921	REAL(wp) :: zrhoy
	922	REAL(wp) :: zcorr
[5656]	923	!!---------------------------------------------
	924	!
[9019]	925	IF( before ) THEN
[636]	926	zrhoy = Agrif_Rhoy()
	927	DO jj=j1,j2
	928	DO ji=i1,i2
[12377]	929	tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj)
[636]	930	END DO
	931	END DO
	932	ELSE
	933	DO jj=j1,j2
	934	DO ji=i1,i2
[9031]	935	tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj)
[4491]	936	!
	937	! Update "now" 3d velocities:
[14064]	938	zpgu(ji,jj) = 0._wp
[4491]	939	DO jk=1,jpkm1
[14064]	940	zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)
[4491]	941	END DO
	942	!
[14064]	943	zcorr = (tabres(ji,jj) - zpgu(ji,jj)) * r1_hu(ji,jj,Kmm_a)
[4491]	944	DO jk=1,jpkm1
[12377]	945	uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)
[4491]	946	END DO
	947	!
[4486]	948	! Update barotropic velocities:
[5930]	949	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
[12489]	950	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
[12377]	951	zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a)
[12489]	952	uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + rn_atfp * zcorr * umask(ji,jj,1)
[5930]	953	END IF
[9031]	954	ENDIF
[12377]	955	uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1)
[4491]	956	!
	957	! Correct "before" velocities to hold correct bt component:
[14064]	958	zpgu(ji,jj) = 0.e0
[4491]	959	DO jk=1,jpkm1
[14064]	960	zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)
[4491]	961	END DO
	962	!
[14064]	963	zcorr = uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a)
[4491]	964	DO jk=1,jpkm1
[12377]	965	uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)
[4491]	966	END DO
	967	!
[636]	968	END DO
	969	END DO
[9031]	970	!
[12489]	971	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	972	uu_b(i1:i2,j1:j2,Kbb_a) = uu_b(i1:i2,j1:j2,Kmm_a)
[9031]	973	ENDIF
[636]	974	ENDIF
[5656]	975	!
[636]	976	END SUBROUTINE updateu2d
[390]	977
[6140]	978
[636]	979	SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
[9019]	980	!!----------------------------------------------------------------------
	981	!! * ROUTINE updatev2d *
	982	!!----------------------------------------------------------------------
	983	INTEGER , INTENT(in ) :: i1, i2, j1, j2
	984	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	985	LOGICAL , INTENT(in ) :: before
[14064]	986	!
	987	REAL(wp), DIMENSION(jpi,jpj) :: zpgv ! 2D workspace
[9019]	988	!
[9031]	989	INTEGER :: ji, jj, jk
[9019]	990	REAL(wp) :: zrhox, zcorr
	991	!!----------------------------------------------------------------------
[5656]	992	!
[9019]	993	IF( before ) THEN
[636]	994	zrhox = Agrif_Rhox()
	995	DO jj=j1,j2
	996	DO ji=i1,i2
[12377]	997	tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj)
[636]	998	END DO
	999	END DO
	1000	ELSE
	1001	DO jj=j1,j2
	1002	DO ji=i1,i2
[9031]	1003	tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj)
[4491]	1004	!
	1005	! Update "now" 3d velocities:
[14064]	1006	zpgv(ji,jj) = 0.e0
[4491]	1007	DO jk=1,jpkm1
[14064]	1008	zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
[4491]	1009	END DO
	1010	!
[14064]	1011	zcorr = (tabres(ji,jj) - zpgv(ji,jj)) * r1_hv(ji,jj,Kmm_a)
[4491]	1012	DO jk=1,jpkm1
[12377]	1013	vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)
[4491]	1014	END DO
	1015	!
[4486]	1016	! Update barotropic velocities:
[5930]	1017	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
[12489]	1018	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
[12377]	1019	zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a)
[12489]	1020	vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + rn_atfp * zcorr * vmask(ji,jj,1)
[5930]	1021	END IF
	1022	ENDIF
[12377]	1023	vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1)
[4491]	1024	!
	1025	! Correct "before" velocities to hold correct bt component:
[14064]	1026	zpgv(ji,jj) = 0.e0
[4491]	1027	DO jk=1,jpkm1
[14064]	1028	zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)
[4491]	1029	END DO
	1030	!
[14064]	1031	zcorr = vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a)
[4491]	1032	DO jk=1,jpkm1
[12377]	1033	vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)
[4491]	1034	END DO
	1035	!
[636]	1036	END DO
	1037	END DO
[9031]	1038	!
[12489]	1039	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	1040	vv_b(i1:i2,j1:j2,Kbb_a) = vv_b(i1:i2,j1:j2,Kmm_a)
[9031]	1041	ENDIF
	1042	!
[636]	1043	ENDIF
[5656]	1044	!
[636]	1045	END SUBROUTINE updatev2d
[390]	1046
[5656]	1047
[636]	1048	SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
[9019]	1049	!!----------------------------------------------------------------------
	1050	!! * ROUTINE updateSSH *
	1051	!!----------------------------------------------------------------------
	1052	INTEGER , INTENT(in ) :: i1, i2, j1, j2
	1053	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	1054	LOGICAL , INTENT(in ) :: before
[5656]	1055	!!
[636]	1056	INTEGER :: ji, jj
[9019]	1057	!!----------------------------------------------------------------------
[5656]	1058	!
[9019]	1059	IF( before ) THEN
[636]	1060	DO jj=j1,j2
	1061	DO ji=i1,i2
[12377]	1062	tabres(ji,jj) = ssh(ji,jj,Kmm_a)
[636]	1063	END DO
	1064	END DO
	1065	ELSE
[12489]	1066	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
[9031]	1067	DO jj=j1,j2
	1068	DO ji=i1,i2
[12377]	1069	ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) &
[12489]	1070	& + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1)
[4491]	1071	END DO
[9031]	1072	END DO
[4491]	1073	ENDIF
[6140]	1074	!
[636]	1075	DO jj=j1,j2
	1076	DO ji=i1,i2
[12377]	1077	ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1)
[636]	1078	END DO
	1079	END DO
[9031]	1080	!
[12489]	1081	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	1082	ssh(i1:i2,j1:j2,Kbb_a) = ssh(i1:i2,j1:j2,Kmm_a)
[9031]	1083	ENDIF
	1084	!
	1085
[636]	1086	ENDIF
[5656]	1087	!
[636]	1088	END SUBROUTINE updateSSH
	1089
[6140]	1090
[4486]	1091	SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
[9019]	1092	!!----------------------------------------------------------------------
	1093	!! * ROUTINE updateub2b *
	1094	!!----------------------------------------------------------------------
	1095	INTEGER , INTENT(in) :: i1, i2, j1, j2
	1096	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	1097	LOGICAL , INTENT(in) :: before
[5656]	1098	!!
[9031]	1099	INTEGER :: ji, jj
	1100	REAL(wp) :: zrhoy, za1, zcor
	1101	!!---------------------------------------------
[5656]	1102	!
[4486]	1103	IF (before) THEN
	1104	zrhoy = Agrif_Rhoy()
	1105	DO jj=j1,j2
	1106	DO ji=i1,i2
	1107	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
	1108	END DO
	1109	END DO
	1110	tabres = zrhoy * tabres
	1111	ELSE
[9031]	1112	!
	1113	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
	1114	!
	1115	za1 = 1._wp / REAL(Agrif_rhot(), wp)
[4486]	1116	DO jj=j1,j2
	1117	DO ji=i1,i2
[9031]	1118	zcor=tabres(ji,jj) - ub2_b(ji,jj)
	1119	! Update time integrated fluxes also in case of multiply nested grids:
	1120	ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor
	1121	! Update corrective fluxes:
	1122	un_bf(ji,jj) = un_bf(ji,jj) + zcor
	1123	! Update half step back fluxes:
	1124	ub2_b(ji,jj) = tabres(ji,jj)
[4486]	1125	END DO
	1126	END DO
	1127	ENDIF
[5656]	1128	!
[4486]	1129	END SUBROUTINE updateub2b
	1130
[9031]	1131	SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
	1132	!!---------------------------------------------
	1133	!! * ROUTINE reflux_sshu *
	1134	!!---------------------------------------------
	1135	INTEGER, INTENT(in) :: i1, i2, j1, j2
	1136	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	1137	LOGICAL, INTENT(in) :: before
	1138	INTEGER, INTENT(in) :: nb, ndir
	1139	!!
	1140	LOGICAL :: western_side, eastern_side
	1141	INTEGER :: ji, jj
	1142	REAL(wp) :: zrhoy, za1, zcor
	1143	!!---------------------------------------------
	1144	!
	1145	IF (before) THEN
	1146	zrhoy = Agrif_Rhoy()
	1147	DO jj=j1,j2
	1148	DO ji=i1,i2
	1149	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
	1150	END DO
	1151	END DO
	1152	tabres = zrhoy * tabres
	1153	ELSE
	1154	!
	1155	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
	1156	!
	1157	western_side = (nb == 1).AND.(ndir == 1)
	1158	eastern_side = (nb == 1).AND.(ndir == 2)
	1159	!
	1160	IF (western_side) THEN
	1161	DO jj=j1,j2
[12489]	1162	zcor = rn_Dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))
[12377]	1163	ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor
[12489]	1164	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + rn_atfp * zcor
[9031]	1165	END DO
	1166	ENDIF
	1167	IF (eastern_side) THEN
	1168	DO jj=j1,j2
[12489]	1169	zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
[12377]	1170	ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor
[12489]	1171	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + rn_atfp * zcor
[9031]	1172	END DO
	1173	ENDIF
	1174	!
	1175	ENDIF
	1176	!
	1177	END SUBROUTINE reflux_sshu
[6140]	1178
[4486]	1179	SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
[9019]	1180	!!----------------------------------------------------------------------
	1181	!! * ROUTINE updatevb2b *
	1182	!!----------------------------------------------------------------------
	1183	INTEGER , INTENT(in ) :: i1, i2, j1, j2
	1184	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	1185	LOGICAL , INTENT(in ) :: before
[5656]	1186	!!
[9031]	1187	INTEGER :: ji, jj
	1188	REAL(wp) :: zrhox, za1, zcor
	1189	!!---------------------------------------------
[5656]	1190	!
[9019]	1191	IF( before ) THEN
[4486]	1192	zrhox = Agrif_Rhox()
	1193	DO jj=j1,j2
	1194	DO ji=i1,i2
	1195	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
	1196	END DO
	1197	END DO
	1198	tabres = zrhox * tabres
	1199	ELSE
[9031]	1200	!
	1201	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
	1202	!
	1203	za1 = 1._wp / REAL(Agrif_rhot(), wp)
[4486]	1204	DO jj=j1,j2
	1205	DO ji=i1,i2
[9031]	1206	zcor=tabres(ji,jj) - vb2_b(ji,jj)
	1207	! Update time integrated fluxes also in case of multiply nested grids:
	1208	vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor
	1209	! Update corrective fluxes:
	1210	vn_bf(ji,jj) = vn_bf(ji,jj) + zcor
	1211	! Update half step back fluxes:
	1212	vb2_b(ji,jj) = tabres(ji,jj)
[4486]	1213	END DO
	1214	END DO
	1215	ENDIF
[5656]	1216	!
[4486]	1217	END SUBROUTINE updatevb2b
	1218
[9031]	1219	SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
	1220	!!---------------------------------------------
	1221	!! * ROUTINE reflux_sshv *
	1222	!!---------------------------------------------
	1223	INTEGER, INTENT(in) :: i1, i2, j1, j2
	1224	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	1225	LOGICAL, INTENT(in) :: before
	1226	INTEGER, INTENT(in) :: nb, ndir
	1227	!!
	1228	LOGICAL :: southern_side, northern_side
	1229	INTEGER :: ji, jj
	1230	REAL(wp) :: zrhox, za1, zcor
	1231	!!---------------------------------------------
	1232	!
	1233	IF (before) THEN
	1234	zrhox = Agrif_Rhox()
	1235	DO jj=j1,j2
	1236	DO ji=i1,i2
	1237	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
	1238	END DO
	1239	END DO
	1240	tabres = zrhox * tabres
	1241	ELSE
	1242	!
	1243	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
	1244	!
	1245	southern_side = (nb == 2).AND.(ndir == 1)
	1246	northern_side = (nb == 2).AND.(ndir == 2)
	1247	!
	1248	IF (southern_side) THEN
	1249	DO ji=i1,i2
[12489]	1250	zcor = rn_Dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1))
[12377]	1251	ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor
[12489]	1252	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + rn_atfp * zcor
[9031]	1253	END DO
	1254	ENDIF
	1255	IF (northern_side) THEN
	1256	DO ji=i1,i2
[12489]	1257	zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2))
[12377]	1258	ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor
[12489]	1259	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + rn_atfp * zcor
[9031]	1260	END DO
	1261	ENDIF
	1262	!
	1263	ENDIF
	1264	!
	1265	END SUBROUTINE reflux_sshv
[5656]	1266
	1267	SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
[6140]	1268	!
[9019]	1269	! ====>>>>>>>>>> currently not used
	1270	!
	1271	!!----------------------------------------------------------------------
	1272	!! * ROUTINE updateT *
	1273	!!----------------------------------------------------------------------
	1274	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	1275	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	1276	LOGICAL , INTENT(in ) :: before
	1277	!!
[5656]	1278	INTEGER :: ji,jj,jk
	1279	REAL(wp) :: ztemp
[9019]	1280	!!----------------------------------------------------------------------
[5656]	1281
	1282	IF (before) THEN
	1283	DO jk=k1,k2
	1284	DO jj=j1,j2
	1285	DO ji=i1,i2
	1286	tabres(ji,jj,jk,1) = e1t(ji,jj)e2t(ji,jj)tmask(ji,jj,jk)
	1287	tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
	1288	tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
	1289	END DO
	1290	END DO
	1291	END DO
	1292	tabres(:,:,:,1)=tabres(:,:,:,1)Agrif_Rhox()Agrif_Rhoy()
	1293	tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
	1294	tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
	1295	ELSE
	1296	DO jk=k1,k2
	1297	DO jj=j1,j2
	1298	DO ji=i1,i2
	1299	IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
	1300	print ,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)e2t(ji,jj)*tmask(ji,jj,jk)
	1301	print ,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)tmask(ji,jj,jk)
	1302	print ,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)tmask(ji,jj,jk)
	1303	ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
	1304	print *,'CORR = ',ztemp-1.
	1305	print ,'NEW VALS = ',tabres(ji,jj,jk,2)ztemp,tabres(ji,jj,jk,3)*ztemp, &
	1306	tabres(ji,jj,jk,2)ztemptabres(ji,jj,jk,3)*ztemp
	1307	e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
	1308	e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
	1309	END IF
	1310	END DO
	1311	END DO
	1312	END DO
	1313	ENDIF
	1314	!
	1315	END SUBROUTINE update_scales
	1316
[6140]	1317
[5656]	1318	SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
[9019]	1319	!!----------------------------------------------------------------------
	1320	!! * ROUTINE updateen *
	1321	!!----------------------------------------------------------------------
	1322	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
	1323	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
	1324	LOGICAL , INTENT(in ) :: before
	1325	!!----------------------------------------------------------------------
[5656]	1326	!
[9019]	1327	IF( before ) THEN
[5656]	1328	ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
	1329	ELSE
	1330	en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
	1331	ENDIF
	1332	!
	1333	END SUBROUTINE updateEN
	1334
	1335
	1336	SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
[9019]	1337	!!----------------------------------------------------------------------
	1338	!! * ROUTINE updateavt *
	1339	!!----------------------------------------------------------------------
	1340	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
	1341	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
	1342	LOGICAL , INTENT(in ) :: before
	1343	!!----------------------------------------------------------------------
[5656]	1344	!
[9019]	1345	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
	1346	ELSE ; avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
[5656]	1347	ENDIF
	1348	!
	1349	END SUBROUTINE updateAVT
	1350
	1351
	1352	SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
	1353	!!---------------------------------------------
	1354	!! * ROUTINE updateavm *
[9019]	1355	!!----------------------------------------------------------------------
	1356	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
	1357	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
	1358	LOGICAL , INTENT(in ) :: before
	1359	!!----------------------------------------------------------------------
[5656]	1360	!
[9019]	1361	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
	1362	ELSE ; avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
[5656]	1363	ENDIF
	1364	!
	1365	END SUBROUTINE updateAVM
	1366
[14053]	1367	#if ! defined key_qco
[9031]	1368	SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
	1369	!!---------------------------------------------
	1370	!! * ROUTINE updatee3t *
	1371	!!---------------------------------------------
	1372	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
	1373	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
	1374	LOGICAL, INTENT(in) :: before
	1375	!
	1376	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
	1377	INTEGER :: ji,jj,jk
	1378	REAL(wp) :: zcoef
	1379	!!---------------------------------------------
	1380	!
	1381	IF (.NOT.before) THEN
	1382	!
	1383	ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))
	1384	!
	1385	! Update e3t from ssh (z* case only)
	1386	DO jk = 1, jpkm1
	1387	DO jj=j1,j2
	1388	DO ji=i1,i2
[12377]	1389	ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) &
[9031]	1390	& *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
	1391	END DO
	1392	END DO
	1393	END DO
	1394	!
	1395	! 1) Updates at BEFORE time step:
	1396	! -------------------------------
	1397	!
	1398	! Save "old" scale factor (prior update) for subsequent asselin correction
	1399	! of prognostic variables
[12377]	1400	e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a)
[9031]	1401
[12377]	1402	! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace...
	1403	! hdiv(i1:i2,j1:j2,1:jpkm1) = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a)
[9031]	1404
[12489]	1405	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN
[9031]	1406	DO jk = 1, jpkm1
	1407	DO jj=j1,j2
	1408	DO ji=i1,i2
[12377]	1409	e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) &
[12489]	1410	& + rn_atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
[9031]	1411	END DO
	1412	END DO
	1413	END DO
	1414	!
[12377]	1415	e3w (i1:i2,j1:j2,1,Kbb_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kbb_a) - e3t_0(i1:i2,j1:j2,1)
	1416	gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp
	1417	gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a)
[9031]	1418	!
	1419	DO jk = 2, jpk
	1420	DO jj = j1,j2
	1421	DO ji = i1,i2
	1422	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
[12377]	1423	e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * &
	1424	& ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) &
[9031]	1425	& + 0.5_wp * tmask(ji,jj,jk) * &
[12377]	1426	& ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) )
	1427	gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a)
	1428	gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) &
	1429	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a))
[9031]	1430	END DO
	1431	END DO
	1432	END DO
	1433	!
	1434	ENDIF
	1435	!
	1436	! 2) Updates at NOW time step:
	1437	! ----------------------------
	1438	!
	1439	! Update vertical scale factor at T-points:
[12377]	1440	e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1)
[9031]	1441	!
	1442	! Update total depth:
[12377]	1443	ht(i1:i2,j1:j2) = 0._wp
[9031]	1444	DO jk = 1, jpkm1
[12377]	1445	ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk)
[9031]	1446	END DO
	1447	!
	1448	! Update vertical scale factor at W-points and depths:
[12377]	1449	e3w (i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1)
	1450	gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a)
	1451	gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp
	1452	gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm_a) - (ht(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
[9031]	1453	!
	1454	DO jk = 2, jpk
	1455	DO jj = j1,j2
	1456	DO ji = i1,i2
	1457	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
[12377]	1458	e3w(ji,jj,jk,Kmm_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) ) &
	1459	& + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) )
	1460	gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a)
	1461	gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) &
	1462	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a))
	1463	gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
[9031]	1464	END DO
	1465	END DO
	1466	END DO
	1467	!
[12489]	1468	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]	1469	e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a)
	1470	e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a)
	1471	gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a)
	1472	gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a)
[9031]	1473	ENDIF
	1474	!
	1475	DEALLOCATE(ptab)
	1476	ENDIF
	1477	!
	1478	END SUBROUTINE updatee3t
[14053]	1479	#endif
[9031]	1480
[390]	1481	#else
[9019]	1482	!!----------------------------------------------------------------------
	1483	!! Empty module no AGRIF zoom
	1484	!!----------------------------------------------------------------------
[636]	1485	CONTAINS
[9570]	1486	SUBROUTINE agrif_oce_update_empty
	1487	WRITE(,) 'agrif_oce_update : You should not have seen this print! error?'
	1488	END SUBROUTINE agrif_oce_update_empty
[390]	1489	#endif
[9019]	1490
	1491	!!======================================================================
[9570]	1492	END MODULE agrif_oce_update
[4486]	1493

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

Download in other formats: