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 NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/sshwzv.F90 @ 14018

Last change on this file since 14018 was 14018, checked in by techene, 3 years ago
#2385 branch updated with trunk 13970
Property svn:keywords set to `Id`
File size: 25.7 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
[13874]	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
	10	!! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work
	11	!! 4.0 ! 2018-12 (A. Coward) add mixed implicit/explicit advection
	12	!! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering.
	13	!! - ! 2020-08 (S. Techene, G. Madec) add here ssh initiatlisation
[3]	14	!!----------------------------------------------------------------------
[1438]	15
[3]	16	!!----------------------------------------------------------------------
[6140]	17	!! ssh_nxt : after ssh
[12377]	18	!! ssh_atf : time filter the ssh arrays
[6140]	19	!! wzv : compute now vertical velocity
[13874]	20	!! ssh_init_rst : ssh set from restart or domcfg.nc file or usr_def_istat_ssh
[1438]	21	!!----------------------------------------------------------------------
[6140]	22	USE oce ! ocean dynamics and tracers variables
[12377]	23	USE isf_oce ! ice shelf
[6140]	24	USE dom_oce ! ocean space and time domain variables
	25	USE sbc_oce ! surface boundary condition: ocean
	26	USE domvvl ! Variable volume
	27	USE divhor ! horizontal divergence
	28	USE phycst ! physical constants
[9019]	29	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
[6140]	30	USE bdydyn2d ! bdy_ssh routine
[2528]	31	#if defined key_agrif
[13286]	32	USE agrif_oce
[9570]	33	USE agrif_oce_interp
[2528]	34	#endif
[6140]	35	!
[10364]	36	USE iom
[6140]	37	USE in_out_manager ! I/O manager
	38	USE restart ! only for lrst_oce
	39	USE prtctl ! Print control
	40	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	41	USE lib_mpp ! MPP library
	42	USE timing ! Timing
[9023]	43	USE wet_dry ! Wetting/Drying flux limiting
[13874]	44	USE usrdef_istate, ONLY : usr_def_istate_ssh ! user defined ssh initial state
	45
[3]	46	IMPLICIT NONE
	47	PRIVATE
	48
[13874]	49	PUBLIC ssh_nxt ! called by step.F90
	50	PUBLIC wzv ! called by step.F90
	51	PUBLIC wAimp ! called by step.F90
	52	PUBLIC ssh_atf ! called by step.F90
	53	PUBLIC ssh_init_rst ! called by domain.F90
[3]	54
	55	!! * Substitutions
[12377]	56	# include "do_loop_substitute.h90"
[13237]	57	# include "domzgr_substitute.h90"
[3]	58	!!----------------------------------------------------------------------
[9598]	59	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	60	!! $Id$
[10068]	61	!! Software governed by the CeCILL license (see ./LICENSE)
[592]	62	!!----------------------------------------------------------------------
[3]	63	CONTAINS
	64
[12377]	65	SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
[3]	66	!!----------------------------------------------------------------------
[4292]	67	!! * ROUTINE ssh_nxt *
[1438]	68	!!
[12377]	69	!! ** Purpose : compute the after ssh (ssh(Kaa))
[3]	70	!!
[4292]	71	!! ** Method : - Using the incompressibility hypothesis, the ssh increment
	72	!! is computed by integrating the horizontal divergence and multiply by
	73	!! by the time step.
[3]	74	!!
[12377]	75	!! ** action : ssh(:,:,Kaa), after sea surface height
[2528]	76	!!
	77	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[3]	78	!!----------------------------------------------------------------------
[12377]	79	INTEGER , INTENT(in ) :: kt ! time step
	80	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
	81	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
[4292]	82	!
[12489]	83	INTEGER :: jk ! dummy loop index
	84	REAL(wp) :: zcoef ! local scalar
[9019]	85	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
[3]	86	!!----------------------------------------------------------------------
[3294]	87	!
[9019]	88	IF( ln_timing ) CALL timing_start('ssh_nxt')
[3294]	89	!
[3]	90	IF( kt == nit000 ) THEN
	91	IF(lwp) WRITE(numout,*)
[4292]	92	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
[1438]	93	IF(lwp) WRITE(numout,*) '~~~~~~~ '
[3]	94	ENDIF
[2528]	95	!
[12489]	96	zcoef = 0.5_wp * r1_rho0
[3]	97
[1438]	98	! !------------------------------!
	99	! ! After Sea Surface Height !
	100	! !------------------------------!
[9023]	101	IF(ln_wd_il) THEN
[12489]	102	CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
[7753]	103	ENDIF
[7646]	104
[12377]	105	CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
[7646]	106	!
[7753]	107	zhdiv(:,:) = 0._wp
[1438]	108	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
[12377]	109	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
[1438]	110	END DO
	111	! ! Sea surface elevation time stepping
[4338]	112	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
	113	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
	114	!
[12489]	115	pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
[9023]	116	!
	117	#if defined key_agrif
[13237]	118	Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa
	119	CALL agrif_ssh( kt )
[9023]	120	#endif
	121	!
[5930]	122	IF ( .NOT.ln_dynspg_ts ) THEN
[7646]	123	IF( ln_bdy ) THEN
[13226]	124	CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
[12377]	125	CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
[5930]	126	ENDIF
[4292]	127	ENDIF
	128	! !------------------------------!
	129	! ! outputs !
	130	! !------------------------------!
	131	!
[12377]	132	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
[4292]	133	!
[9019]	134	IF( ln_timing ) CALL timing_stop('ssh_nxt')
[4292]	135	!
	136	END SUBROUTINE ssh_nxt
	137
	138
[13237]	139	SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
[4292]	140	!!----------------------------------------------------------------------
	141	!! * ROUTINE wzv *
	142	!!
	143	!! ** Purpose : compute the now vertical velocity
	144	!!
	145	!! ** Method : - Using the incompressibility hypothesis, the vertical
	146	!! velocity is computed by integrating the horizontal divergence
	147	!! from the bottom to the surface minus the scale factor evolution.
	148	!! The boundary conditions are w=0 at the bottom (no flux) and.
	149	!!
[12377]	150	!! ** action : pww : now vertical velocity
[4292]	151	!!
	152	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	153	!!----------------------------------------------------------------------
[12377]	154	INTEGER , INTENT(in) :: kt ! time step
	155	INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices
[13237]	156	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
[4292]	157	!
[5836]	158	INTEGER :: ji, jj, jk ! dummy loop indices
[9019]	159	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
[4292]	160	!!----------------------------------------------------------------------
	161	!
[9019]	162	IF( ln_timing ) CALL timing_start('wzv')
[5836]	163	!
[4292]	164	IF( kt == nit000 ) THEN
	165	IF(lwp) WRITE(numout,*)
	166	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
	167	IF(lwp) WRITE(numout,*) '~~~~~ '
	168	!
[12377]	169	pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
[4292]	170	ENDIF
	171	! !------------------------------!
	172	! ! Now Vertical Velocity !
	173	! !------------------------------!
	174	!
[13237]	175	! !===============================!
	176	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
	177	! !===============================!
[9019]	178	ALLOCATE( zhdiv(jpi,jpj,jpk) )
[4292]	179	!
	180	DO jk = 1, jpkm1
	181	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
[4338]	182	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
[13295]	183	DO_2D( 0, 0, 0, 0 )
[12377]	184	zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
	185	END_2D
[592]	186	END DO
[13226]	187	CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
[4292]	188	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
[12377]	189	! ! Same question holds for hdiv. Perhaps just for security
[4292]	190	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	191	! computation of w
[13237]	192	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
	193	& + zhdiv(:,:,jk) &
	194	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
	195	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
[4292]	196	END DO
[12377]	197	! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
[9019]	198	DEALLOCATE( zhdiv )
[13237]	199	! !=================================!
	200	ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
	201	! !=================================!
	202	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	203	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) ) * tmask(:,:,jk)
	204	END DO
	205	! !==========================================!
	206	ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
	207	! !==========================================!
[13998]	208	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
[12377]	209	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
[13237]	210	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
	211	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
[4292]	212	END DO
[1438]	213	ENDIF
[592]	214
[7646]	215	IF( ln_bdy ) THEN
[4327]	216	DO jk = 1, jpkm1
[12377]	217	pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
[4327]	218	END DO
	219	ENDIF
[4292]	220	!
[13286]	221	#if defined key_agrif
	222	IF( .NOT. AGRIF_Root() ) THEN
	223	!
[12965]	224	! Mask vertical velocity at first/last columns/row
	225	! inside computational domain (cosmetic)
[13286]	226	DO jk = 1, jpkm1
	227	IF( lk_west ) THEN ! --- West --- !
	228	DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
	229	DO jj = 1, jpj
	230	pww(ji,jj,jk) = 0._wp
	231	END DO
	232	END DO
	233	ENDIF
	234	IF( lk_east ) THEN ! --- East --- !
	235	DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
	236	DO jj = 1, jpj
	237	pww(ji,jj,jk) = 0._wp
	238	END DO
	239	END DO
	240	ENDIF
	241	IF( lk_south ) THEN ! --- South --- !
	242	DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
	243	DO ji = 1, jpi
	244	pww(ji,jj,jk) = 0._wp
	245	END DO
	246	END DO
	247	ENDIF
	248	IF( lk_north ) THEN ! --- North --- !
	249	DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
	250	DO ji = 1, jpi
	251	pww(ji,jj,jk) = 0._wp
	252	END DO
	253	END DO
	254	ENDIF
	255	!
	256	END DO
[12965]	257	!
[9023]	258	ENDIF
[13286]	259	#endif
[5836]	260	!
[9124]	261	IF( ln_timing ) CALL timing_stop('wzv')
[9023]	262	!
[5836]	263	END SUBROUTINE wzv
[592]	264
	265
[13237]	266	SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh, pssh_f )
[1438]	267	!!----------------------------------------------------------------------
[12377]	268	!! * ROUTINE ssh_atf *
[1438]	269	!!
[12377]	270	!! ** Purpose : Apply Asselin time filter to now SSH.
[1438]	271	!!
[2528]	272	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
	273	!! from the filter, see Leclair and Madec 2010) and swap :
[12489]	274	!! pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
	275	!! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
[1438]	276	!!
[12377]	277	!! ** action : - pssh(:,:,Kmm) time filtered
[2528]	278	!!
	279	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]	280	!!----------------------------------------------------------------------
[12377]	281	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	282	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
[13237]	283	REAL(wp), DIMENSION(jpi,jpj,jpt) , TARGET, INTENT(inout) :: pssh ! SSH field
	284	REAL(wp), DIMENSION(jpi,jpj ), OPTIONAL, TARGET, INTENT( out) :: pssh_f ! filtered SSH field
[6140]	285	!
	286	REAL(wp) :: zcoef ! local scalar
[13237]	287	REAL(wp), POINTER, DIMENSION(:,:) :: zssh ! pointer for filtered SSH
[1438]	288	!!----------------------------------------------------------------------
[3294]	289	!
[12377]	290	IF( ln_timing ) CALL timing_start('ssh_atf')
[3294]	291	!
[1438]	292	IF( kt == nit000 ) THEN
	293	IF(lwp) WRITE(numout,*)
[12377]	294	IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
[1438]	295	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	296	ENDIF
[6140]	297	! !== Euler time-stepping: no filter, just swap ==!
[12489]	298	IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps
[13237]	299	IF( PRESENT( pssh_f ) ) THEN ; zssh => pssh_f
	300	ELSE ; zssh => pssh(:,:,Kmm)
	301	ENDIF
[12377]	302	! ! filtered "now" field
[12489]	303	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
[13427]	304	!
[12377]	305	IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed
[12489]	306	zcoef = rn_atfp * rn_Dt * r1_rho0
[12377]	307	pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) &
	308	& - rnf_b(:,:) + rnf (:,:) &
	309	& + fwfisf_cav_b(:,:) - fwfisf_cav(:,:) &
	310	& + fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * ssmask(:,:)
	311
	312	! ice sheet coupling
[13427]	313	IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 ) &
	314	& pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
[12377]	315
[6140]	316	ENDIF
[1438]	317	ENDIF
	318	!
[13427]	319	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf - pssh(:,:,Kmm): ', mask1=tmask )
[2528]	320	!
[12377]	321	IF( ln_timing ) CALL timing_stop('ssh_atf')
[3294]	322	!
[12377]	323	END SUBROUTINE ssh_atf
[3]	324
[13237]	325
[12377]	326	SUBROUTINE wAimp( kt, Kmm )
[10364]	327	!!----------------------------------------------------------------------
	328	!! * ROUTINE wAimp *
	329	!!
	330	!! ** Purpose : compute the Courant number and partition vertical velocity
	331	!! if a proportion needs to be treated implicitly
	332	!!
	333	!! ** Method : -
	334	!!
[12377]	335	!! ** action : ww : now vertical velocity (to be handled explicitly)
[10364]	336	!! : wi : now vertical velocity (for implicit treatment)
	337	!!
[11414]	338	!! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
	339	!! implicit scheme for vertical advection in oceanic modeling.
	340	!! Ocean Modelling, 91, 38-69.
[10364]	341	!!----------------------------------------------------------------------
	342	INTEGER, INTENT(in) :: kt ! time step
[12377]	343	INTEGER, INTENT(in) :: Kmm ! time level index
[10364]	344	!
	345	INTEGER :: ji, jj, jk ! dummy loop indices
[13237]	346	REAL(wp) :: zCu, zcff, z1_e3t, zdt ! local scalars
[10364]	347	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
[11407]	348	REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters
[10364]	349	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
	350	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
	351	!!----------------------------------------------------------------------
	352	!
	353	IF( ln_timing ) CALL timing_start('wAimp')
	354	!
	355	IF( kt == nit000 ) THEN
	356	IF(lwp) WRITE(numout,*)
	357	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
	358	IF(lwp) WRITE(numout,*) '~~~~~ '
[11293]	359	wi(:,:,:) = 0._wp
[10364]	360	ENDIF
	361	!
[11414]	362	! Calculate Courant numbers
[13237]	363	zdt = 2._wp * rn_Dt ! 2*rn_Dt and not rDt (for restartability)
[11414]	364	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
[13295]	365	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	366	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[13237]	367	Cu_adv(ji,jj,jk) = zdt * &
	368	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
	369	& + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) &
	370	& * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - &
	371	& MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) &
	372	& * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) &
[12377]	373	& * r1_e1e2t(ji,jj) &
[13237]	374	& + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) &
	375	& * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - &
	376	& MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) &
	377	& * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) &
[12377]	378	& * r1_e1e2t(ji,jj) &
	379	& ) * z1_e3t
	380	END_3D
[11414]	381	ELSE
[13295]	382	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	383	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[13237]	384	Cu_adv(ji,jj,jk) = zdt * &
	385	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
[12377]	386	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm), 0._wp ) - &
	387	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm), 0._wp ) ) &
	388	& * r1_e1e2t(ji,jj) &
	389	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm), 0._wp ) - &
	390	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm), 0._wp ) ) &
	391	& * r1_e1e2t(ji,jj) &
	392	& ) * z1_e3t
	393	END_3D
[11414]	394	ENDIF
[13226]	395	CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
[10364]	396	!
	397	CALL iom_put("Courant",Cu_adv)
	398	!
	399	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
[13998]	400	DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! or scan Courant criterion and partition ! w where necessary
[12377]	401	!
	402	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
[11293]	403	! alt:
[12377]	404	! IF ( ww(ji,jj,jk) > 0._wp ) THEN
[11293]	405	! zCu = Cu_adv(ji,jj,jk)
	406	! ELSE
	407	! zCu = Cu_adv(ji,jj,jk-1)
	408	! ENDIF
[12377]	409	!
	410	IF( zCu <= Cu_min ) THEN !<-- Fully explicit
	411	zcff = 0._wp
	412	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
	413	zcff = ( zCu - Cu_min )**2
	414	zcff = zcff / ( Fcu + zcff )
	415	ELSE !<-- Mostly implicit
	416	zcff = ( zCu - Cu_max )/ zCu
	417	ENDIF
	418	zcff = MIN(1._wp, zcff)
	419	!
	420	wi(ji,jj,jk) = zcff * ww(ji,jj,jk)
	421	ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
	422	!
	423	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl
	424	END_3D
[11293]	425	Cu_adv(:,:,1) = 0._wp
[10364]	426	ELSE
	427	! Fully explicit everywhere
[11407]	428	Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl
[10907]	429	wi (:,:,:) = 0._wp
[10364]	430	ENDIF
	431	CALL iom_put("wimp",wi)
	432	CALL iom_put("wi_cff",Cu_adv)
[12377]	433	CALL iom_put("wexp",ww)
[10364]	434	!
	435	IF( ln_timing ) CALL timing_stop('wAimp')
	436	!
	437	END SUBROUTINE wAimp
[13874]	438
	439
	440	SUBROUTINE ssh_init_rst( Kbb, Kmm, Kaa )
	441	!!---------------------------------------------------------------------
	442	!! * ROUTINE ssh_init_rst *
	443	!!
	444	!! ** Purpose : ssh initialization of the sea surface height (ssh)
	445	!!
	446	!! ** Method : set ssh from restart or read configuration, or user_def
	447	!! * ln_rstart = T
	448	!! USE of IOM library to read ssh in the restart file
	449	!! Leap-Frog: Kbb and Kmm are read except for l_1st_euler=T
[13895]	450	!!
[13874]	451	!! * otherwise
	452	!! call user defined ssh or
	453	!! set to -ssh_ref in wet and drying case with domcfg.nc
	454	!!
	455	!! NB: ssh_b/n are written by restart.F90
	456	!!----------------------------------------------------------------------
	457	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
	458	!
	459	INTEGER :: ji, jj, jk
	460	!!----------------------------------------------------------------------
	461	!
[13915]	462	IF(lwp) THEN
	463	WRITE(numout,*)
	464	WRITE(numout,*) 'ssh_init_rst : ssh initialization'
	465	WRITE(numout,*) '~~~~~~~~~~~~ '
	466	ENDIF
	467	!
[13874]	468	! !=============================!
	469	IF( ln_rstart ) THEN !== Read the restart file ==!
	470	! !=============================!
	471	!
[13895]	472	! !* Read ssh at Kmm
[13915]	473	IF(lwp) WRITE(numout,*)
	474	IF(lwp) WRITE(numout,*) ' Kmm sea surface height read in the restart file'
[14018]	475	CALL iom_get( numror, jpdom_auto, 'sshn' , ssh(:,:,Kmm) )
[13874]	476	!
[13895]	477	IF( l_1st_euler ) THEN !* Euler at first time-step
[13915]	478	IF(lwp) WRITE(numout,*)
	479	IF(lwp) WRITE(numout,*) ' Euler first time step : ssh(Kbb) = ssh(Kmm)'
[13895]	480	ssh(:,:,Kbb) = ssh(:,:,Kmm)
[13874]	481	!
[13895]	482	ELSE !* read ssh at Kbb
[13915]	483	IF(lwp) WRITE(numout,*)
	484	IF(lwp) WRITE(numout,*) ' Kbb sea surface height read in the restart file'
[14018]	485	CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb) )
[13874]	486	ENDIF
	487	! !============================!
	488	ELSE !== Initialize at "rest" ==!
	489	! !============================!
	490	!
[13915]	491	IF(lwp) WRITE(numout,*)
	492	IF(lwp) WRITE(numout,*) ' initialization at rest'
	493	!
[13874]	494	IF( ll_wd ) THEN !* wet and dry
	495	!
	496	IF( ln_read_cfg ) THEN ! read configuration : ssh_ref is read in domain_cfg file
	497	!!st why ssh is not masked : i.e. ssh(:,:,Kmm) = -ssh_ref*ssmask(:,:),
	498	!!st since at the 1st time step lbclnk will be applied on ssh at Kaa but not initially at Kbb and Kmm
	499	ssh(:,:,Kbb) = -ssh_ref
	500	!
	501	DO_2D( 1, 1, 1, 1 )
	502	IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth
	503	ssh(ji,jj,Kbb) = rn_wdmin1 - ht_0(ji,jj)
	504	ENDIF
	505	END_2D
	506	ELSE ! user define configuration case
	507	CALL usr_def_istate_ssh( tmask, ssh(:,:,Kbb) )
	508	ENDIF
	509	!
	510	ELSE !* user defined configuration
	511	CALL usr_def_istate_ssh( tmask, ssh(:,:,Kbb) )
	512	!
	513	ENDIF
	514	!
	515	ssh(:,:,Kmm) = ssh(:,:,Kbb) !* set now values from to before ones
	516	ssh(:,:,Kaa) = 0._wp
	517	ENDIF
	518	!
	519	END SUBROUTINE ssh_init_rst
	520
[3]	521	!!======================================================================
[1565]	522	END MODULE sshwzv

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

Download in other formats: