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.

sshwzv.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 4427

Last change on this file since 4427 was 4427, checked in by trackstand2, 10 years ago
First files changed on last FINISS work package. Stephen's work although commited by Andy P.
Property svn:keywords set to `Id`
File size: 25.5 KB

Rev	Line
[1565]	1	MODULE sshwzv
[3]	2	!!==============================================================================
[1438]	3	!! * MODULE sshwzv *
	4	!! Ocean dynamics : sea surface height and vertical velocity
[3]	5	!!==============================================================================
[1438]	6	!! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code
[2528]	7	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA
	8	!! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
	9	!! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module
[3]	10	!!----------------------------------------------------------------------
[1438]	11
[3]	12	!!----------------------------------------------------------------------
[1438]	13	!! ssh_wzv : after ssh & now vertical velocity
	14	!! ssh_nxt : filter ans swap the ssh arrays
	15	!!----------------------------------------------------------------------
[3]	16	USE oce ! ocean dynamics and tracers variables
	17	USE dom_oce ! ocean space and time domain variables
[888]	18	USE sbc_oce ! surface boundary condition: ocean
	19	USE domvvl ! Variable volume
[1565]	20	USE divcur ! hor. divergence and curl (div & cur routines)
[1438]	21	USE iom ! I/O library
	22	USE restart ! only for lrst_oce
[3]	23	USE in_out_manager ! I/O manager
[258]	24	USE prtctl ! Print control
[592]	25	USE phycst
	26	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[2715]	27	USE lib_mpp ! MPP library
[1241]	28	USE obc_par ! open boundary cond. parameter
	29	USE obc_oce
[2528]	30	USE bdy_oce
[2715]	31	USE diaar5, ONLY: lk_diaar5
[1482]	32	USE iom
[2715]	33	USE sbcrnf, ONLY: h_rnf, nk_rnf ! River runoff
[2528]	34	#if defined key_agrif
	35	USE agrif_opa_update
[2486]	36	USE agrif_opa_interp
[2528]	37	#endif
	38	#if defined key_asminc
	39	USE asminc ! Assimilation increment
	40	#endif
[592]	41
[3]	42	IMPLICIT NONE
	43	PRIVATE
	44
[1438]	45	PUBLIC ssh_wzv ! called by step.F90
	46	PUBLIC ssh_nxt ! called by step.F90
[3]	47
[3211]	48	!! * Control permutation of array indices
	49	# include "oce_ftrans.h90"
	50	# include "dom_oce_ftrans.h90"
	51	# include "sbc_oce_ftrans.h90"
	52	# include "domvvl_ftrans.h90"
	53	# include "obc_oce_ftrans.h90"
	54	#if defined key_asminc
	55	# include "asminc_ftrans.h90"
	56	#endif
	57
[3]	58	!! * Substitutions
	59	# include "domzgr_substitute.h90"
[1438]	60	# include "vectopt_loop_substitute.h90"
[3]	61	!!----------------------------------------------------------------------
[2528]	62	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	63	!! $Id$
[2715]	64	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[592]	65	!!----------------------------------------------------------------------
[3]	66	CONTAINS
	67
[1438]	68	SUBROUTINE ssh_wzv( kt )
[3]	69	!!----------------------------------------------------------------------
[1438]	70	!! * ROUTINE ssh_wzv *
	71	!!
	72	!! ** Purpose : compute the after ssh (ssha), the now vertical velocity
	73	!! and update the now vertical coordinate (lk_vvl=T).
[3]	74	!!
[2528]	75	!! ** Method : - Using the incompressibility hypothesis, the vertical
[1438]	76	!! velocity is computed by integrating the horizontal divergence
	77	!! from the bottom to the surface minus the scale factor evolution.
	78	!! The boundary conditions are w=0 at the bottom (no flux) and.
[3]	79	!!
[1438]	80	!! ** action : ssha : after sea surface height
	81	!! wn : now vertical velocity
[2528]	82	!! sshu_a, sshv_a, sshf_a : after sea surface height (lk_vvl=T)
	83	!! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
	84	!!
	85	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[3]	86	!!----------------------------------------------------------------------
[2715]	87	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	88	USE oce , ONLY: z3d => ta ! ta used as 3D workspace
	89	USE wrk_nemo, ONLY: zhdiv => wrk_2d_1 , z2d => wrk_2d_2 ! 2D workspace
[3211]	90	!! DCSE_NEMO: need additional directives for renamed module variables
	91	!FTRANS z3d :I :I :z
[2715]	92	!
[1438]	93	INTEGER, INTENT(in) :: kt ! time step
[2715]	94	!
	95	INTEGER :: ji, jj, jk ! dummy loop indices
[4427]	96	#if defined key_z_first
	97	INTEGER :: klim ! upper bound on k loop
	98	#endif
[2715]	99	REAL(wp) :: zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0 ! local scalars
[3]	100	!!----------------------------------------------------------------------
	101
[2715]	102	IF( wrk_in_use(2, 1,2) ) THEN
	103	CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable') ; RETURN
	104	ENDIF
	105
[3]	106	IF( kt == nit000 ) THEN
[2528]	107	!
[3]	108	IF(lwp) WRITE(numout,*)
[1438]	109	IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
	110	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	111	!
[4427]	112	#if defined key_z_first
	113	DO jj=1,jpj
	114	DO ji=1,jpi
	115	DO jk=mbkmax(ji,jj), jpk
	116	wn(ji,jj,jk) = 0._wp
	117	END DO
	118	END DO
	119	END DO
	120	#else
[2715]	121	wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
[4427]	122	#endif
[1438]	123	!
	124	IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only)
	125	DO jj = 1, jpjm1
	126	DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible
[3211]	127	#if defined key_z_first
	128	zcoefu = 0.5 * umask_1(ji,jj) / ( e1u(ji,jj) * e2u(ji,jj) )
	129	zcoefv = 0.5 * vmask_1(ji,jj) / ( e1v(ji,jj) * e2v(ji,jj) )
	130	zcoeff = 0.25 * umask_1(ji,jj) * umask_1(ji,jj+1)
	131	#else
[1438]	132	zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
	133	zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
	134	zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
[3211]	135	#endif
[1438]	136	sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &
	137	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
	138	sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &
	139	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
	140	sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &
	141	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
	142	sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &
	143	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
	144	END DO
	145	END DO
	146	CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )
	147	CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )
[2528]	148	DO jj = 1, jpjm1
	149	DO ji = 1, jpim1 ! NO Vector Opt.
[3211]	150	#if defined key_z_first
	151	sshf_n(ji,jj) = 0.5 * umask_1(ji,jj) * umask_1(ji,jj+1) &
	152	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
	153	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
	154	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
	155	#else
[2528]	156	sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &
	157	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
	158	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
	159	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
[3211]	160	#endif
[2528]	161	END DO
	162	END DO
	163	CALL lbc_lnk( sshf_n, 'F', 1. )
[1438]	164	ENDIF
	165	!
[3]	166	ENDIF
	167
[2528]	168	! !------------------------------------------!
	169	IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case)
	170	! !------------------------------------------!
[3211]	171	#if defined key_z_first
[3432]	172	! DCSE_NEMO: can't use implicit loop over k here because the domzgr_substitute.h90
	173	! file causes the line below to be expanded to:
	174	! gdept_1(1:jpkm1,:,:) = (gdept(1:jpkm1,:,:)(1.+sshn(:,:)mut(1:jpkm1,:,:)))
	175	! which contains non-conforming array expressions.
[4427]	176	DO jj=1,jpj
	177	DO ji=1,jpi
	178	klim=mbkmax(ji,jj)
	179	! now local depths stored in fsdep. arrays
	180	fsdept(ji,jj,1:klim) = fsdept_n(ji,jj,1:klim)
	181	fsdepw(ji,jj,1:klim) = fsdepw_n(ji,jj,1:klim)
	182	fsde3w(ji,jj,1:klim) = fsde3w_n(ji,jj,1:klim)
	183	! vertical scale factors stored in fse3. arrays
	184	fse3t (ji,jj,1:klim) = fse3t_n (ji,jj,1:klim)
	185	fse3u (ji,jj,1:klim) = fse3u_n (ji,jj,1:klim)
	186	fse3v (ji,jj,1:klim) = fse3v_n (ji,jj,1:klim)
	187	fse3f (ji,jj,1:klim) = fse3f_n (ji,jj,1:klim)
	188	fse3w (ji,jj,1:klim) = fse3w_n (ji,jj,1:klim)
	189	fse3uw(ji,jj,1:klim) = fse3uw_n(ji,jj,1:klim)
	190	fse3vw(ji,jj,1:klim) = fse3vw_n(ji,jj,1:klim)
[3432]	191	END DO
	192	END DO
[3211]	193	#else
[1565]	194	DO jk = 1, jpkm1
[2528]	195	fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays
[1565]	196	fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
	197	fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
	198	!
[2528]	199	fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays
[1565]	200	fse3u (:,:,jk) = fse3u_n (:,:,jk)
	201	fse3v (:,:,jk) = fse3v_n (:,:,jk)
	202	fse3f (:,:,jk) = fse3f_n (:,:,jk)
	203	fse3w (:,:,jk) = fse3w_n (:,:,jk)
	204	fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
	205	fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
	206	END DO
[3211]	207	#endif
[2528]	208	!
	209	hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points)
[1565]	210	hv(:,:) = hv_0(:,:) + sshv_n(:,:)
[2528]	211	! ! now masked inverse of the ocean depth (at u- and v-points)
[3211]	212	#if defined key_z_first
	213	hur(:,:) = umask_1(:,:) / ( hu(:,:) + 1._wp - umask_1(:,:) )
	214	hvr(:,:) = vmask_1(:,:) / ( hv(:,:) + 1._wp - vmask_1(:,:) )
	215	#else
[2715]	216	hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
	217	hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
[3211]	218	#endif
[2528]	219	!
[1565]	220	ENDIF
[2528]	221	!
	222	CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity
	223	!
[2715]	224	z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog)
	225	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
[3]	226
[4415]	227	WRITE(,)'ARPDBG, ssh_wzv: sum WWW of hdivn=',SUM(hdivn),' at step=',kt
	228	WRITE(,)'ARPDBG, ssh_wzv: sum WWW of fse3t=',SUM(fse3t(:,:,:)),' at step=',kt
	229
[1438]	230	! !------------------------------!
	231	! ! After Sea Surface Height !
	232	! !------------------------------!
[2715]	233	zhdiv(:,:) = 0._wp
[3211]	234	#if defined key_z_first
	235	DO jj = 1, jpj
	236	DO ji = 1, jpi
[4427]	237	DO jk = 1, mbkmax(ji,jj)-1 ! Horizontal divergence of barotropic transports
[3211]	238	zhdiv(ji,jj) = zhdiv(ji,jj) + fse3t(ji,jj,jk) * hdivn(ji,jj,jk)
	239	END DO
	240	END DO
	241	END DO
	242	#else
[1438]	243	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
	244	zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
	245	END DO
[3211]	246	#endif
[4415]	247	WRITE(,)'ARPDBG, ssh_wzv: sum XXX of zhdiv=',SUM(zhdiv),' at step=',kt
	248	WRITE(,)'ARPDBG, ssh_wzv: sum XXX of ssha=',SUM(ssha),' at step=',kt
[1438]	249	! ! Sea surface elevation time stepping
[2528]	250	! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
	251	! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
	252	z1_rau0 = 0.5 / rau0
[3211]	253	#if defined key_z_first
	254	ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask_1(:,:)
	255	#else
[2715]	256	ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1)
[3211]	257	#endif
[4415]	258	WRITE(,)'ARPDBG, ssh_wzv: sum YYY of sshb=',SUM(sshb),' at step=',kt
	259	WRITE(,)'ARPDBG, ssh_wzv: sum YYY of ssha=',SUM(ssha),' at step=',kt
[2486]	260	#if defined key_agrif
[2715]	261	CALL agrif_ssh( kt )
[2486]	262	#endif
[1438]	263	#if defined key_obc
[2528]	264	IF( Agrif_Root() ) THEN
[1438]	265	ssha(:,:) = ssha(:,:) * obctmsk(:,:)
[2715]	266	CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm)
[1438]	267	ENDIF
	268	#endif
[2528]	269	#if defined key_bdy
	270	ssha(:,:) = ssha(:,:) * bdytmask(:,:)
	271	CALL lbc_lnk( ssha, 'T', 1. )
	272	#endif
[1438]	273
	274	! ! Sea Surface Height at u-,v- and f-points (vvl case only)
	275	IF( lk_vvl ) THEN ! (required only in key_vvl case)
	276	DO jj = 1, jpjm1
[1694]	277	DO ji = 1, jpim1 ! NO Vector Opt.
[3211]	278	#if defined key_z_first
	279	sshu_a(ji,jj) = 0.5 * umask_1(ji,jj) / ( e1u(ji ,jj) * e2u(ji ,jj) ) &
	280	& * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) &
	281	& + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
	282	sshv_a(ji,jj) = 0.5 * vmask_1(ji,jj) / ( e1v(ji,jj ) * e2v(ji,jj ) ) &
	283	& * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) &
	284	& + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
	285	#else
[1438]	286	sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) &
	287	& * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) &
	288	& + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
	289	sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) &
	290	& * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) &
	291	& + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
[3211]	292	#endif
[592]	293	END DO
	294	END DO
[2715]	295	CALL lbc_lnk( sshu_a, 'U', 1. ) ; CALL lbc_lnk( sshv_a, 'V', 1. ) ! Boundaries conditions
[1438]	296	ENDIF
[2715]	297
[2528]	298	#if defined key_asminc
[2715]	299	! ! Include the IAU weighted SSH increment
	300	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
[2528]	301	CALL ssh_asm_inc( kt )
	302	ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
	303	ENDIF
	304	#endif
[592]	305
[1438]	306	! !------------------------------!
	307	! ! Now Vertical Velocity !
	308	! !------------------------------!
[2528]	309	z1_2dt = 1.e0 / z2dt
[3211]	310	#if defined key_z_first
	311	DO jj = 1, jpj
	312	DO ji = 1, jpi
[4427]	313	DO jk = mbkmax(ji,jj)-1, 1, -1 ! integrate from the bottom the hor. divergence
[3211]	314	wn(ji,jj,jk) = wn(ji,jj,jk+1) &
	315	& - fse3t_n(ji,jj,jk) * hdivn(ji,jj,jk) &
	316	& - ( fse3t_a(ji,jj,jk) - fse3t_b(ji,jj,jk) ) &
	317	& * tmask(ji,jj,jk) * z1_2dt
	318	#if defined key_bdy
	319	wn(ji,jj,jk) = wn(ji,jj,jk) * bdytmask(ji,jj)
	320	#endif
	321	END DO
	322	END DO
	323	END DO
	324	#else
[2528]	325	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	326	! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
[2715]	327	wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) &
	328	& - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) &
[2528]	329	& * tmask(:,:,jk) * z1_2dt
	330	#if defined key_bdy
	331	wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
	332	#endif
[1438]	333	END DO
[3211]	334	#endif
[2528]	335
	336	! !------------------------------!
	337	! ! outputs !
	338	! !------------------------------!
[1756]	339	CALL iom_put( "woce", wn ) ! vertical velocity
	340	CALL iom_put( "ssh" , sshn ) ! sea surface height
	341	CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height
[2528]	342	IF( lk_diaar5 ) THEN ! vertical mass transport & its square value
	343	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[1756]	344	z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
[3211]	345	#if defined key_z_first
	346	DO jj = 1, jpj
	347	DO ji = 1, jpi
[4427]	348	DO jk = 1, mbkmax(ji,jj)
[3211]	349	z3d(ji,jj,jk) = wn(ji,jj,jk) * z2d(ji,jj)
	350	END DO
	351	END DO
	352	END DO
	353	#else
[1756]	354	DO jk = 1, jpk
	355	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
	356	END DO
[3211]	357	#endif
[2528]	358	CALL iom_put( "w_masstr" , z3d )
	359	CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
[1756]	360	ENDIF
[1438]	361	!
[2528]	362	IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 )
	363	!
[2715]	364	IF( wrk_not_released(2, 1,2) ) CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
	365	!
[1438]	366	END SUBROUTINE ssh_wzv
[592]	367
	368
[1438]	369	SUBROUTINE ssh_nxt( kt )
	370	!!----------------------------------------------------------------------
	371	!! * ROUTINE ssh_nxt *
	372	!!
	373	!! ** Purpose : achieve the sea surface height time stepping by
	374	!! applying Asselin time filter and swapping the arrays
	375	!! ssha already computed in ssh_wzv
	376	!!
[2528]	377	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
	378	!! from the filter, see Leclair and Madec 2010) and swap :
	379	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
	380	!! - atfp * rdt * ( emp_b - emp ) / rau0
	381	!! sshn = ssha
[1438]	382	!!
	383	!! ** action : - sshb, sshn : before & now sea surface height
	384	!! ready for the next time step
[2528]	385	!!
	386	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]	387	!!----------------------------------------------------------------------
[2528]	388	INTEGER, INTENT(in) :: kt ! ocean time-step index
[1438]	389	!!
[2528]	390	INTEGER :: ji, jj ! dummy loop indices
	391	REAL(wp) :: zec ! temporary scalar
[1438]	392	!!----------------------------------------------------------------------
[592]	393
[1438]	394	IF( kt == nit000 ) THEN
	395	IF(lwp) WRITE(numout,*)
	396	IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
	397	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	398	ENDIF
[592]	399
[2528]	400	! !--------------------------!
[2715]	401	IF( lk_vvl ) THEN ! Variable volume levels ! (ssh at t-, u-, v, f-points)
[2528]	402	! !--------------------------!
	403	!
[2715]	404	IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter
	405	sshn (:,:) = ssha (:,:) ! now <-- after (before already = now)
[1438]	406	sshu_n(:,:) = sshu_a(:,:)
	407	sshv_n(:,:) = sshv_a(:,:)
[2715]	408	DO jj = 1, jpjm1 ! ssh now at f-point
[2528]	409	DO ji = 1, jpim1 ! NO Vector Opt.
[3211]	410	#if defined key_z_first
	411	sshf_n(ji,jj) = 0.5 * umask_1(ji,jj) * umask_1(ji,jj+1) &
	412	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
	413	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
	414	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
	415	#else
[2715]	416	sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &
[2528]	417	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
	418	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
	419	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
[3211]	420	#endif
[2528]	421	END DO
	422	END DO
[2715]	423	CALL lbc_lnk( sshf_n, 'F', 1. ) ! Boundaries conditions
	424	!
	425	ELSE !** Leap-Frog time-stepping: Asselin filter + swap
[4415]	426	WRITE(,) 'ARPDBG: ssh_nxt: SUM of sshb = ',SUM(sshb),' at step=',kt
	427	WRITE(,) 'ARPDBG: ssh_nxt: SUM of sshn = ',SUM(sshn),' at step=',kt
[2528]	428	zec = atfp * rdt / rau0
[1438]	429	DO jj = 1, jpj
[2715]	430	DO ji = 1, jpi ! before <-- now filtered
[3211]	431	#if defined key_z_first
[2528]	432	sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) &
[3211]	433	& - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask_1(ji,jj)
	434	#else
	435	sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) &
[2528]	436	& - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
[3211]	437	#endif
[2715]	438	sshn (ji,jj) = ssha (ji,jj) ! now <-- after
[1438]	439	sshu_n(ji,jj) = sshu_a(ji,jj)
	440	sshv_n(ji,jj) = sshv_a(ji,jj)
	441	END DO
	442	END DO
[2715]	443	DO jj = 1, jpjm1 ! ssh now at f-point
[2528]	444	DO ji = 1, jpim1 ! NO Vector Opt.
[3211]	445	#if defined key_z_first
	446	sshf_n(ji,jj) = 0.5 * umask_1(ji,jj) * umask_1(ji,jj+1) &
	447	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
	448	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
	449	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
	450	#else
[2528]	451	sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &
	452	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
	453	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
	454	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
[3211]	455	#endif
[2528]	456	END DO
	457	END DO
[2715]	458	CALL lbc_lnk( sshf_n, 'F', 1. ) ! Boundaries conditions
	459	!
	460	DO jj = 1, jpjm1 ! ssh before at u- & v-points
[2528]	461	DO ji = 1, jpim1 ! NO Vector Opt.
[3211]	462	#if defined key_z_first
	463	sshu_b(ji,jj) = 0.5 * umask_1(ji,jj) / ( e1u(ji ,jj) * e2u(ji ,jj) ) &
	464	& * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &
	465	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
	466	sshv_b(ji,jj) = 0.5 * vmask_1(ji,jj) / ( e1v(ji,jj ) * e2v(ji,jj ) ) &
	467	& * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &
	468	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
	469	#else
[2528]	470	sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) &
	471	& * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &
	472	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
	473	sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) &
	474	& * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &
	475	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
[3211]	476	#endif
[2528]	477	END DO
	478	END DO
	479	CALL lbc_lnk( sshu_b, 'U', 1. )
[2715]	480	CALL lbc_lnk( sshv_b, 'V', 1. ) ! Boundaries conditions
	481	!
[1438]	482	ENDIF
[2528]	483	! !--------------------------!
[2715]	484	ELSE ! fixed levels ! (ssh at t-point only)
[2528]	485	! !--------------------------!
[1438]	486	!
[2715]	487	IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter
	488	sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)
[1438]	489	!
[2715]	490	ELSE ! Leap-Frog time-stepping: Asselin filter + swap
[1438]	491	DO jj = 1, jpj
[2715]	492	DO ji = 1, jpi ! before <-- now filtered
[2528]	493	sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
[2715]	494	sshn(ji,jj) = ssha(ji,jj) ! now <-- after
[1438]	495	END DO
	496	END DO
	497	ENDIF
	498	!
	499	ENDIF
	500	!
[2528]	501	! Update velocity at AGRIF zoom boundaries
[2486]	502	#if defined key_agrif
[2528]	503	IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
[2486]	504	#endif
[1438]	505	!
[2528]	506	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )
	507	!
[1438]	508	END SUBROUTINE ssh_nxt
[3]	509
	510	!!======================================================================
[1565]	511	END MODULE sshwzv

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

Download in other formats: