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/devukmo2010/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/devukmo2010/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2128

Last change on this file since 2128 was 2128, checked in by rfurner, 14 years ago
merged branches OBS, ASM, Rivers, BDY & mixed_dynldf ready for vn3.3 merge
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 15.6 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
[2128]	7	!! 3.3 ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
	8	!! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module
[3]	9	!!----------------------------------------------------------------------
[1438]	10
[3]	11	!!----------------------------------------------------------------------
[1438]	12	!! ssh_wzv : after ssh & now vertical velocity
	13	!! ssh_nxt : filter ans swap the ssh arrays
	14	!!----------------------------------------------------------------------
[3]	15	USE oce ! ocean dynamics and tracers variables
	16	USE dom_oce ! ocean space and time domain variables
[888]	17	USE sbc_oce ! surface boundary condition: ocean
	18	USE domvvl ! Variable volume
[1565]	19	USE divcur ! hor. divergence and curl (div & cur routines)
	20	USE cla_div ! cross land: hor. divergence (div_cla routine)
[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)
[1241]	27	USE obc_par ! open boundary cond. parameter
	28	USE obc_oce
[2128]	29	USE bdy_oce
[1756]	30	USE diaar5, ONLY : lk_diaar5
[1482]	31	USE iom
[2128]	32	USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep ! River runoff
	33	#if defined key_asminc
	34	USE asminc ! Assimilation increment
	35	#endif
[592]	36
[3]	37	IMPLICIT NONE
	38	PRIVATE
	39
[1438]	40	PUBLIC ssh_wzv ! called by step.F90
	41	PUBLIC ssh_nxt ! called by step.F90
[3]	42
	43	!! * Substitutions
	44	# include "domzgr_substitute.h90"
[1438]	45	# include "vectopt_loop_substitute.h90"
	46
[3]	47	!!----------------------------------------------------------------------
[1438]	48	!! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
[888]	49	!! $Id$
[592]	50	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
	51	!!----------------------------------------------------------------------
[3]	52
	53	CONTAINS
	54
[1438]	55	SUBROUTINE ssh_wzv( kt )
[3]	56	!!----------------------------------------------------------------------
[1438]	57	!! * ROUTINE ssh_wzv *
	58	!!
	59	!! ** Purpose : compute the after ssh (ssha), the now vertical velocity
	60	!! and update the now vertical coordinate (lk_vvl=T).
[3]	61	!!
[1438]	62	!! ** Method : -
[3]	63	!!
[1438]	64	!! - Using the incompressibility hypothesis, the vertical
	65	!! velocity is computed by integrating the horizontal divergence
	66	!! from the bottom to the surface minus the scale factor evolution.
	67	!! The boundary conditions are w=0 at the bottom (no flux) and.
[3]	68	!!
[1438]	69	!! ** action : ssha : after sea surface height
	70	!! wn : now vertical velocity
	71	!! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height
	72	!! at u-, v-, f-point s
	73	!! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
[3]	74	!!----------------------------------------------------------------------
[1756]	75	USE oce, ONLY : z3d => ta ! use ta as 3D workspace
	76	!!
[1438]	77	INTEGER, INTENT(in) :: kt ! time step
	78	!!
	79	INTEGER :: ji, jj, jk ! dummy loop indices
	80	REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars
	81	REAL(wp) :: z2dt, zraur ! temporary scalars
	82	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
[1756]	83	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
[3]	84	!!----------------------------------------------------------------------
	85
	86	IF( kt == nit000 ) THEN
	87	IF(lwp) WRITE(numout,*)
[1438]	88	IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
	89	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	90	!
	91	wn(:,:,jpk) = 0.e0 ! bottom boundary condition: w=0 (set once for all)
	92	!
	93	IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only)
	94	DO jj = 1, jpjm1
	95	DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible
	96	zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
	97	zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
	98	zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
	99	sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &
	100	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
	101	sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &
	102	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
	103	sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) &
	104	& + sshb(ji+1,jj) + sshb(ji+1,jj+1) )
	105	sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &
	106	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
	107	sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &
	108	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
	109	sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) &
	110	& + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
	111	END DO
	112	END DO
	113	CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )
	114	CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )
	115	CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. )
	116	ENDIF
	117	!
[3]	118	ENDIF
	119
[1565]	120	! !------------------------------!
[1607]	121	IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case)
[1565]	122	! !------------------------------!
	123	DO jk = 1, jpkm1
[1607]	124	fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays
[1565]	125	fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
	126	fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
	127	!
[1607]	128	fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays
[1565]	129	fse3u (:,:,jk) = fse3u_n (:,:,jk)
	130	fse3v (:,:,jk) = fse3v_n (:,:,jk)
	131	fse3f (:,:,jk) = fse3f_n (:,:,jk)
	132	fse3w (:,:,jk) = fse3w_n (:,:,jk)
	133	fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
	134	fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
	135	END DO
[1607]	136	! ! now ocean depth (at u- and v-points)
[1565]	137	hu(:,:) = hu_0(:,:) + sshu_n(:,:)
	138	hv(:,:) = hv_0(:,:) + sshv_n(:,:)
[1607]	139	! ! now masked inverse of the ocean depth (at u- and v-points)
[1565]	140	hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) )
	141	hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
	142	!
[2128]	143	DO jj=1,jpj
	144	DO ji=1,jpi
	145	rnf_dep(ji,jj)=0
	146	DO jk=1,rnf_mod_dep(ji,jj) ! recalculates rnf_dep to be the depth
	147	rnf_dep(ji,jj)=rnf_dep(ji,jj)+fse3t(ji,jj,jk) ! in metres to the bottom of the relevant grid box
	148	ENDDO
	149	ENDDO
	150	ENDDO
	151	!
[1565]	152	ENDIF
	153
	154	CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity
	155	IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence)
	156
[1438]	157	! set time step size (Euler/Leapfrog)
	158	z2dt = 2. * rdt
[1607]	159	IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt
[3]	160
[1739]	161	zraur = 1. / rau0
[592]	162
[1438]	163	! !------------------------------!
	164	! ! After Sea Surface Height !
	165	! !------------------------------!
	166	zhdiv(:,:) = 0.e0
	167	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
	168	zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
	169	END DO
	170
	171	! ! Sea surface elevation time stepping
	172	ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)
	173
	174	#if defined key_obc
	175	IF ( Agrif_Root() ) THEN
	176	ssha(:,:) = ssha(:,:) * obctmsk(:,:)
	177	CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm)
	178	ENDIF
	179	#endif
	180
[2128]	181	#if defined key_bdy
	182	ssha(:,:) = ssha(:,:) * bdytmask(:,:)
	183	CALL lbc_lnk( ssha, 'T', 1. )
	184	#endif
	185
[1438]	186	! ! Sea Surface Height at u-,v- and f-points (vvl case only)
	187	IF( lk_vvl ) THEN ! (required only in key_vvl case)
	188	DO jj = 1, jpjm1
[1694]	189	DO ji = 1, jpim1 ! NO Vector Opt.
[1438]	190	sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) &
	191	& * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) &
	192	& + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
	193	sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) &
	194	& * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) &
	195	& + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
[1607]	196	sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) &
[1438]	197	& * ( ssha(ji ,jj) + ssha(ji ,jj+1) &
	198	& + ssha(ji+1,jj) + ssha(ji+1,jj+1) )
[592]	199	END DO
	200	END DO
[1438]	201	CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions
	202	CALL lbc_lnk( sshv_a, 'V', 1. )
	203	CALL lbc_lnk( sshf_a, 'F', 1. )
	204	ENDIF
[592]	205
[2128]	206	! Include the IAU weighted SSH increment
	207	#if defined key_asminc
	208	IF( ( lk_asminc ).AND.( ln_sshinc ).AND.( ln_asmiau ) ) THEN
	209	CALL ssh_asm_inc( kt )
	210	ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
	211	ENDIF
	212	#endif
	213
[1438]	214	! !------------------------------!
	215	! ! Now Vertical Velocity !
	216	! !------------------------------!
	217	! ! integrate from the bottom the hor. divergence
	218	DO jk = jpkm1, 1, -1
	219	wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) &
	220	& - ( fse3t_a(:,:,jk) &
	221	& - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt
[2128]	222	#if defined key_bdy
	223	wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
	224	#endif
[1438]	225	END DO
[1482]	226	!
[1756]	227	CALL iom_put( "woce", wn ) ! vertical velocity
	228	CALL iom_put( "ssh" , sshn ) ! sea surface height
	229	CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height
	230	IF( lk_diaar5 ) THEN
	231	z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
	232	DO jk = 1, jpk
	233	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
	234	END DO
	235	CALL iom_put( "w_masstr" , z3d ) ! vertical mass transport
	236	CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) ! square of vertical mass transport
	237	ENDIF
[1438]	238	!
	239	END SUBROUTINE ssh_wzv
[592]	240
	241
[1438]	242	SUBROUTINE ssh_nxt( kt )
	243	!!----------------------------------------------------------------------
	244	!! * ROUTINE ssh_nxt *
	245	!!
	246	!! ** Purpose : achieve the sea surface height time stepping by
	247	!! applying Asselin time filter and swapping the arrays
	248	!! ssha already computed in ssh_wzv
	249	!!
	250	!! ** Method : - apply Asselin time fiter to now ssh and swap :
	251	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
	252	!! sshn = ssha
	253	!!
	254	!! ** action : - sshb, sshn : before & now sea surface height
	255	!! ready for the next time step
	256	!!----------------------------------------------------------------------
	257	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	258	!!
	259	INTEGER :: ji, jj ! dummy loop indices
	260	!!----------------------------------------------------------------------
[592]	261
[1438]	262	IF( kt == nit000 ) THEN
	263	IF(lwp) WRITE(numout,*)
	264	IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
	265	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	266	ENDIF
[592]	267
[1438]	268	! Time filter and swap of the ssh
	269	! -------------------------------
[592]	270
[1438]	271	IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points
	272	! ! ---------------------- !
	273	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter
	274	sshn (:,:) = ssha (:,:) ! now <-- after (before already = now)
	275	sshu_n(:,:) = sshu_a(:,:)
	276	sshv_n(:,:) = sshv_a(:,:)
	277	sshf_n(:,:) = sshf_a(:,:)
	278	ELSE ! Leap-Frog time-stepping: Asselin filter + swap
	279	DO jj = 1, jpj
	280	DO ji = 1, jpi ! before <-- now filtered
	281	sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) )
	282	sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) )
	283	sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) )
	284	sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) )
	285	sshn (ji,jj) = ssha (ji,jj) ! now <-- after
	286	sshu_n(ji,jj) = sshu_a(ji,jj)
	287	sshv_n(ji,jj) = sshv_a(ji,jj)
	288	sshf_n(ji,jj) = sshf_a(ji,jj)
	289	END DO
	290	END DO
	291	ENDIF
	292	!
	293	ELSE ! fixed levels : ssh at t-point only
	294	! ! ------------ !
	295	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter
	296	sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)
	297	!
	298	ELSE ! Leap-Frog time-stepping: Asselin filter + swap
	299	DO jj = 1, jpj
	300	DO ji = 1, jpi ! before <-- now filtered
	301	sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
	302	sshn(ji,jj) = ssha(ji,jj) ! now <-- after
	303	END DO
	304	END DO
	305	ENDIF
	306	!
	307	ENDIF
	308	!
	309	IF(ln_ctl) CALL prt_ctl(tab2d_1=sshb , clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )
	310	!
	311	END SUBROUTINE ssh_nxt
[3]	312
	313	!!======================================================================
[1565]	314	END MODULE sshwzv

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

Download in other formats: