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/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 5221

Last change on this file since 5221 was 5014, checked in by hliu, 10 years ago
upload the modifications for W/D based on r:4826
Property svn:keywords set to `Id`
File size: 14.2 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
[4292]	10	!! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work
[5014]	11	!! 3.? ! 2014-09 (H. Liu) add wetting and drying
[3]	12	!!----------------------------------------------------------------------
[1438]	13
[3]	14	!!----------------------------------------------------------------------
[4292]	15	!! ssh_nxt : after ssh
	16	!! ssh_swp : filter ans swap the ssh arrays
	17	!! wzv : compute now vertical velocity
[1438]	18	!!----------------------------------------------------------------------
[3]	19	USE oce ! ocean dynamics and tracers variables
	20	USE dom_oce ! ocean space and time domain variables
[888]	21	USE sbc_oce ! surface boundary condition: ocean
	22	USE domvvl ! Variable volume
[1565]	23	USE divcur ! hor. divergence and curl (div & cur routines)
[1438]	24	USE iom ! I/O library
[4292]	25	USE restart ! only for lrst_oce
[3]	26	USE in_out_manager ! I/O manager
[258]	27	USE prtctl ! Print control
[592]	28	USE phycst
	29	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[2715]	30	USE lib_mpp ! MPP library
[2528]	31	USE bdy_oce
[4292]	32	USE bdy_par
	33	USE bdydyn2d ! bdy_ssh routine
[2715]	34	USE diaar5, ONLY: lk_diaar5
[1482]	35	USE iom
[2528]	36	#if defined key_agrif
	37	USE agrif_opa_update
[2486]	38	USE agrif_opa_interp
[2528]	39	#endif
	40	#if defined key_asminc
	41	USE asminc ! Assimilation increment
	42	#endif
[3294]	43	USE wrk_nemo ! Memory Allocation
	44	USE timing ! Timing
[5014]	45	USE wadlmt ! Wetting/Drying flux limting
[592]	46
[3]	47	IMPLICIT NONE
	48	PRIVATE
	49
[1438]	50	PUBLIC ssh_nxt ! called by step.F90
[4292]	51	PUBLIC wzv ! called by step.F90
	52	PUBLIC ssh_swp ! called by step.F90
[3]	53
	54	!! * Substitutions
	55	# include "domzgr_substitute.h90"
[1438]	56	# include "vectopt_loop_substitute.h90"
[3]	57	!!----------------------------------------------------------------------
[2528]	58	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	59	!! $Id$
[2715]	60	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[592]	61	!!----------------------------------------------------------------------
[3]	62	CONTAINS
	63
[4292]	64	SUBROUTINE ssh_nxt( kt )
[3]	65	!!----------------------------------------------------------------------
[4292]	66	!! * ROUTINE ssh_nxt *
[1438]	67	!!
[4292]	68	!! ** Purpose : compute the after ssh (ssha)
[3]	69	!!
[4292]	70	!! ** Method : - Using the incompressibility hypothesis, the ssh increment
	71	!! is computed by integrating the horizontal divergence and multiply by
	72	!! by the time step.
[3]	73	!!
[1438]	74	!! ** action : ssha : after sea surface height
[2528]	75	!!
	76	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[3]	77	!!----------------------------------------------------------------------
[2715]	78	!
[4292]	79	REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv
	80	INTEGER, INTENT(in) :: kt ! time step
	81	!
	82	INTEGER :: jk ! dummy loop indice
	83	REAL(wp) :: z2dt, z1_rau0 ! local scalars
[3]	84	!!----------------------------------------------------------------------
[3294]	85	!
[4292]	86	IF( nn_timing == 1 ) CALL timing_start('ssh_nxt')
[3294]	87	!
[4292]	88	CALL wrk_alloc( jpi, jpj, zhdiv )
[3294]	89	!
[3]	90	IF( kt == nit000 ) THEN
[2528]	91	!
[3]	92	IF(lwp) WRITE(numout,*)
[4292]	93	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
[1438]	94	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	95	!
[3]	96	ENDIF
[2528]	97	!
[2715]	98	z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog)
	99	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
[5014]	100	!
	101	z1_rau0 = 0.5_wp * r1_rau0
	102	!
	103	IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt)
[3]	104
[5014]	105	CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity
	106	!
	107
[1438]	108	! !------------------------------!
	109	! ! After Sea Surface Height !
	110	! !------------------------------!
[2715]	111	zhdiv(:,:) = 0._wp
[1438]	112	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
[4292]	113	zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
[1438]	114	END DO
	115	! ! Sea surface elevation time stepping
[4338]	116	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
	117	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
	118	!
[2715]	119	ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1)
[1438]	120
[4486]	121	#if ! defined key_dynspg_ts
	122	! These lines are not necessary with time splitting since
	123	! boundary condition on sea level is set during ts loop
[2486]	124	#if defined key_agrif
[2715]	125	CALL agrif_ssh( kt )
[2486]	126	#endif
[2528]	127	#if defined key_bdy
[4292]	128	IF (lk_bdy) THEN
	129	CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
	130	CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
	131	ENDIF
[2528]	132	#endif
[4486]	133	#endif
	134
[3764]	135	#if defined key_asminc
	136	! ! Include the IAU weighted SSH increment
	137	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
	138	CALL ssh_asm_inc( kt )
	139	ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
	140	ENDIF
	141	#endif
[4292]	142
	143	! !------------------------------!
	144	! ! outputs !
	145	! !------------------------------!
	146	CALL iom_put( "ssh" , sshn ) ! sea surface height
	147	CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height
	148	!
	149	IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 )
	150	!
	151	CALL wrk_dealloc( jpi, jpj, zhdiv )
	152	!
	153	IF( nn_timing == 1 ) CALL timing_stop('ssh_nxt')
	154	!
	155	END SUBROUTINE ssh_nxt
	156
	157
	158	SUBROUTINE wzv( kt )
	159	!!----------------------------------------------------------------------
	160	!! * ROUTINE wzv *
	161	!!
	162	!! ** Purpose : compute the now vertical velocity
	163	!!
	164	!! ** Method : - Using the incompressibility hypothesis, the vertical
	165	!! velocity is computed by integrating the horizontal divergence
	166	!! from the bottom to the surface minus the scale factor evolution.
	167	!! The boundary conditions are w=0 at the bottom (no flux) and.
	168	!!
	169	!! ** action : wn : now vertical velocity
	170	!!
	171	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	172	!!----------------------------------------------------------------------
	173	!
	174	INTEGER, INTENT(in) :: kt ! time step
	175	REAL(wp), POINTER, DIMENSION(:,: ) :: z2d
	176	REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv
	177	!
	178	INTEGER :: ji, jj, jk ! dummy loop indices
	179	REAL(wp) :: z1_2dt ! local scalars
	180	!!----------------------------------------------------------------------
	181
	182	IF( nn_timing == 1 ) CALL timing_start('wzv')
	183	!
	184	IF( kt == nit000 ) THEN
	185	!
	186	IF(lwp) WRITE(numout,*)
	187	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
	188	IF(lwp) WRITE(numout,*) '~~~~~ '
	189	!
	190	wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
	191	!
	192	ENDIF
	193	! !------------------------------!
	194	! ! Now Vertical Velocity !
	195	! !------------------------------!
	196	z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog)
	197	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt
	198	!
	199	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases
	200	CALL wrk_alloc( jpi, jpj, jpk, zhdiv )
	201	!
	202	DO jk = 1, jpkm1
	203	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
[4338]	204	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
[4292]	205	DO jj = 2, jpjm1
	206	DO ji = fs_2, fs_jpim1 ! vector opt.
	207	zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
	208	END DO
[592]	209	END DO
	210	END DO
[4292]	211	CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
	212	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
	213	! ! Same question holds for hdivn. Perhaps just for security
	214	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	215	! computation of w
	216	wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) &
	217	& + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
	218	END DO
	219	! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
	220	CALL wrk_dealloc( jpi, jpj, jpk, zhdiv )
	221	ELSE ! z_star and linear free surface cases
	222	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	223	! computation of w
	224	wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) &
	225	& + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
	226	END DO
[1438]	227	ENDIF
[592]	228
[2528]	229	#if defined key_bdy
[4327]	230	IF (lk_bdy) THEN
	231	DO jk = 1, jpkm1
	232	wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
	233	END DO
	234	ENDIF
[2528]	235	#endif
[4292]	236	!
[2528]	237	! !------------------------------!
	238	! ! outputs !
	239	! !------------------------------!
[4292]	240	CALL iom_put( "woce", wn ) ! vertical velocity
[2528]	241	IF( lk_diaar5 ) THEN ! vertical mass transport & its square value
[4292]	242	CALL wrk_alloc( jpi, jpj, z2d )
	243	CALL wrk_alloc( jpi, jpj, jpk, z3d )
[2528]	244	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[4292]	245	z2d(:,:) = rau0 * e12t(:,:)
[1756]	246	DO jk = 1, jpk
	247	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
	248	END DO
[2528]	249	CALL iom_put( "w_masstr" , z3d )
	250	CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
[4292]	251	CALL wrk_dealloc( jpi, jpj, z2d )
	252	CALL wrk_dealloc( jpi, jpj, jpk, z3d )
[1756]	253	ENDIF
[1438]	254	!
[4292]	255	IF( nn_timing == 1 ) CALL timing_stop('wzv')
[592]	256
	257
[4292]	258	END SUBROUTINE wzv
	259
	260	SUBROUTINE ssh_swp( kt )
[1438]	261	!!----------------------------------------------------------------------
	262	!! * ROUTINE ssh_nxt *
	263	!!
	264	!! ** Purpose : achieve the sea surface height time stepping by
	265	!! applying Asselin time filter and swapping the arrays
[4292]	266	!! ssha already computed in ssh_nxt
[1438]	267	!!
[2528]	268	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
	269	!! from the filter, see Leclair and Madec 2010) and swap :
	270	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
	271	!! - atfp * rdt * ( emp_b - emp ) / rau0
	272	!! sshn = ssha
[1438]	273	!!
	274	!! ** action : - sshb, sshn : before & now sea surface height
	275	!! ready for the next time step
[2528]	276	!!
	277	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]	278	!!----------------------------------------------------------------------
[2528]	279	INTEGER, INTENT(in) :: kt ! ocean time-step index
[1438]	280	!!----------------------------------------------------------------------
[3294]	281	!
[4292]	282	IF( nn_timing == 1 ) CALL timing_start('ssh_swp')
[3294]	283	!
[1438]	284	IF( kt == nit000 ) THEN
	285	IF(lwp) WRITE(numout,*)
[4292]	286	IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
[1438]	287	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	288	ENDIF
[592]	289
[4292]	290	# if defined key_dynspg_ts
	291	IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter
	292	# else
	293	IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter
	294	#endif
	295	sshb(:,:) = sshn(:,:) ! before <-- now
	296	sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)
	297	ELSE !** Leap-Frog time-stepping: Asselin filter + swap
	298	sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered
	299	IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
	300	sshn(:,:) = ssha(:,:) ! now <-- after
[1438]	301	ENDIF
	302	!
[2528]	303	! Update velocity at AGRIF zoom boundaries
[2486]	304	#if defined key_agrif
[2528]	305	IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
[2486]	306	#endif
[1438]	307	!
[2528]	308	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )
	309	!
[4292]	310	IF( nn_timing == 1 ) CALL timing_stop('ssh_swp')
[3294]	311	!
[4292]	312	END SUBROUTINE ssh_swp
[3]	313
	314	!!======================================================================
[1565]	315	END MODULE sshwzv

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

Download in other formats: