Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sshwzv.F90 @ 11738

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 11738

Download in other formats: