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_r12377_KERNEL-06_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90 @ 12680

Last change on this file since 12680 was 12680, checked in by techene, 4 years ago
dynatfQCO.F90, stepLF.F90 : fixed (remove pe3. from dyn_atf_qco input arguments), all : remove e3. tables and include gurvan's feedbacks
Property svn:keywords set to `Id`
File size: 19.9 KB

Rev	Line
[12590]	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
[12590]	7	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA
[2528]	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
[11414]	11	!! 4.0 ! 2018-12 (A. Coward) add mixed implicit/explicit advection
[12377]	12	!! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering.
[3]	13	!!----------------------------------------------------------------------
[1438]	14
[3]	15	!!----------------------------------------------------------------------
[6140]	16	!! ssh_nxt : after ssh
[12377]	17	!! ssh_atf : time filter the ssh arrays
[6140]	18	!! wzv : compute now vertical velocity
[1438]	19	!!----------------------------------------------------------------------
[6140]	20	USE oce ! ocean dynamics and tracers variables
[12377]	21	USE isf_oce ! ice shelf
[12590]	22	USE dom_oce ! ocean space and time domain variables
[6140]	23	USE sbc_oce ! surface boundary condition: ocean
	24	USE domvvl ! Variable volume
	25	USE divhor ! horizontal divergence
	26	USE phycst ! physical constants
[9019]	27	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
[6140]	28	USE bdydyn2d ! bdy_ssh routine
[2528]	29	#if defined key_agrif
[9570]	30	USE agrif_oce_interp
[2528]	31	#endif
[6140]	32	!
[12590]	33	USE iom
[6140]	34	USE in_out_manager ! I/O manager
	35	USE restart ! only for lrst_oce
	36	USE prtctl ! Print control
	37	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	38	USE lib_mpp ! MPP library
	39	USE timing ! Timing
[9023]	40	USE wet_dry ! Wetting/Drying flux limiting
[592]	41
[3]	42	IMPLICIT NONE
	43	PRIVATE
	44
[1438]	45	PUBLIC ssh_nxt ! called by step.F90
[4292]	46	PUBLIC wzv ! called by step.F90
[10364]	47	PUBLIC wAimp ! called by step.F90
[12377]	48	PUBLIC ssh_atf ! called by step.F90
[3]	49
	50	!! * Substitutions
[12377]	51	# include "do_loop_substitute.h90"
[12590]	52	# include "domzgr_substitute.h90"
	53
[3]	54	!!----------------------------------------------------------------------
[9598]	55	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	56	!! $Id$
[10068]	57	!! Software governed by the CeCILL license (see ./LICENSE)
[592]	58	!!----------------------------------------------------------------------
[3]	59	CONTAINS
	60
[12377]	61	SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
[3]	62	!!----------------------------------------------------------------------
[4292]	63	!! * ROUTINE ssh_nxt *
[12590]	64	!!
[12377]	65	!! ** Purpose : compute the after ssh (ssh(Kaa))
[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	!!
[12377]	71	!! ** action : ssh(:,:,Kaa), after sea surface height
[2528]	72	!!
	73	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[3]	74	!!----------------------------------------------------------------------
[12377]	75	INTEGER , INTENT(in ) :: kt ! time step
	76	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
	77	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
[12590]	78	!
[7753]	79	INTEGER :: jk ! dummy loop indice
[5836]	80	REAL(wp) :: z2dt, zcoef ! local scalars
[9019]	81	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
[3]	82	!!----------------------------------------------------------------------
[3294]	83	!
[9019]	84	IF( ln_timing ) CALL timing_start('ssh_nxt')
[3294]	85	!
[3]	86	IF( kt == nit000 ) THEN
	87	IF(lwp) WRITE(numout,*)
[4292]	88	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
[1438]	89	IF(lwp) WRITE(numout,*) '~~~~~~~ '
[3]	90	ENDIF
[2528]	91	!
[7646]	92	z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog)
[2715]	93	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
[7646]	94	zcoef = 0.5_wp * r1_rau0
[3]	95
[1438]	96	! !------------------------------!
	97	! ! After Sea Surface Height !
	98	! !------------------------------!
[9023]	99	IF(ln_wd_il) THEN
[12377]	100	CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv )
[7753]	101	ENDIF
[7646]	102
[12377]	103	CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
[7646]	104	!
[7753]	105	zhdiv(:,:) = 0._wp
[1438]	106	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
[12377]	107	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
[1438]	108	END DO
	109	! ! Sea surface elevation time stepping
[4338]	110	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
	111	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
[12590]	112	!
[12377]	113	pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
[9023]	114	!
	115	#if defined key_agrif
[12680]	116	Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa
	117	CALL agrif_ssh( kt )
[9023]	118	#endif
	119	!
[5930]	120	IF ( .NOT.ln_dynspg_ts ) THEN
[7646]	121	IF( ln_bdy ) THEN
[12377]	122	CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. ) ! Not sure that's necessary
	123	CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
[5930]	124	ENDIF
[4292]	125	ENDIF
	126	! !------------------------------!
	127	! ! outputs !
	128	! !------------------------------!
	129	!
[12377]	130	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
[4292]	131	!
[9019]	132	IF( ln_timing ) CALL timing_stop('ssh_nxt')
[4292]	133	!
	134	END SUBROUTINE ssh_nxt
	135
[12590]	136
[12680]	137	SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
[4292]	138	!!----------------------------------------------------------------------
	139	!! * ROUTINE wzv *
[12590]	140	!!
[4292]	141	!! ** Purpose : compute the now vertical velocity
	142	!!
[12590]	143	!! ** Method : - Using the incompressibility hypothesis, the vertical
	144	!! velocity is computed by integrating the horizontal divergence
[4292]	145	!! from the bottom to the surface minus the scale factor evolution.
	146	!! The boundary conditions are w=0 at the bottom (no flux) and.
	147	!!
[12377]	148	!! ** action : pww : now vertical velocity
[4292]	149	!!
	150	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	151	!!----------------------------------------------------------------------
[12377]	152	INTEGER , INTENT(in) :: kt ! time step
	153	INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices
[12680]	154	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
[4292]	155	!
[5836]	156	INTEGER :: ji, jj, jk ! dummy loop indices
	157	REAL(wp) :: z1_2dt ! local scalars
[9019]	158	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
[4292]	159	!!----------------------------------------------------------------------
	160	!
[9019]	161	IF( ln_timing ) CALL timing_start('wzv')
[5836]	162	!
[4292]	163	IF( kt == nit000 ) THEN
	164	IF(lwp) WRITE(numout,*)
	165	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
	166	IF(lwp) WRITE(numout,*) '~~~~~ '
	167	!
[12680]	168	pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
[4292]	169	ENDIF
[12680]	170	!
	171	z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog)
[4292]	172	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt
	173	!
[12680]	174	! !===============================!
	175	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
	176	! !===============================!
[12590]	177	ALLOCATE( zhdiv(jpi,jpj,jpk) )
[4292]	178	!
	179	DO jk = 1, jpkm1
	180	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
[4338]	181	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
[12377]	182	DO_2D_00_00
	183	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) )
	184	END_2D
[592]	185	END DO
[10425]	186	CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
[4292]	187	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
[12377]	188	! ! Same question holds for hdiv. Perhaps just for security
[4292]	189	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	190	! computation of w
[12680]	191	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
	192	& + zhdiv(:,:,jk) &
	193	& + z1_2dt * ( e3t(:,:,jk,Kaa) &
	194	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
[4292]	195	END DO
[12377]	196	! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
[12590]	197	DEALLOCATE( zhdiv )
[12680]	198	! !=================================!
	199	ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
	200	! !=================================!
	201	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
	202	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) ) * tmask(:,:,jk)
	203	END DO
	204	! !==========================================!
	205	ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
	206	! !==========================================!
	207	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
[12590]	208	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
	209	& + z1_2dt * ( e3t(:,:,jk,Kaa) &
	210	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
[4292]	211	END DO
[1438]	212	ENDIF
[592]	213
[7646]	214	IF( ln_bdy ) THEN
[4327]	215	DO jk = 1, jpkm1
[12377]	216	pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
[4327]	217	END DO
	218	ENDIF
[4292]	219	!
[12590]	220	#if defined key_agrif
	221	IF( .NOT. AGRIF_Root() ) THEN
	222	IF ((nbondi == 1).OR.(nbondi == 2)) pww(nlci-1 , : ,:) = 0.e0 ! east
	223	IF ((nbondi == -1).OR.(nbondi == 2)) pww(2 , : ,:) = 0.e0 ! west
	224	IF ((nbondj == 1).OR.(nbondj == 2)) pww(: ,nlcj-1 ,:) = 0.e0 ! north
	225	IF ((nbondj == -1).OR.(nbondj == 2)) pww(: ,2 ,:) = 0.e0 ! south
	226	ENDIF
	227	#endif
[5836]	228	!
[9124]	229	IF( ln_timing ) CALL timing_stop('wzv')
[9023]	230	!
[5836]	231	END SUBROUTINE wzv
[592]	232
	233
[12377]	234	SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
[1438]	235	!!----------------------------------------------------------------------
[12377]	236	!! * ROUTINE ssh_atf *
[1438]	237	!!
[12377]	238	!! ** Purpose : Apply Asselin time filter to now SSH.
[1438]	239	!!
[2528]	240	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
	241	!! from the filter, see Leclair and Madec 2010) and swap :
[12377]	242	!! pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
[2528]	243	!! - atfp * rdt * ( emp_b - emp ) / rau0
[1438]	244	!!
[12377]	245	!! ** action : - pssh(:,:,Kmm) time filtered
[2528]	246	!!
	247	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]	248	!!----------------------------------------------------------------------
[12377]	249	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	250	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
	251	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field
[6140]	252	!
	253	REAL(wp) :: zcoef ! local scalar
[1438]	254	!!----------------------------------------------------------------------
[3294]	255	!
[12377]	256	IF( ln_timing ) CALL timing_start('ssh_atf')
[3294]	257	!
[1438]	258	IF( kt == nit000 ) THEN
	259	IF(lwp) WRITE(numout,*)
[12377]	260	IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
[1438]	261	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	262	ENDIF
[6140]	263	! !== Euler time-stepping: no filter, just swap ==!
[12377]	264	IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN ! Only do time filtering for leapfrog timesteps
	265	! ! filtered "now" field
	266	pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
	267	IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed
[6140]	268	zcoef = atfp * rdt * r1_rau0
[12377]	269	pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) &
	270	& - rnf_b(:,:) + rnf (:,:) &
	271	& + fwfisf_cav_b(:,:) - fwfisf_cav(:,:) &
	272	& + fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * ssmask(:,:)
	273
	274	! ice sheet coupling
	275	IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
	276
[6140]	277	ENDIF
[1438]	278	ENDIF
	279	!
[12377]	280	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm) - : ', mask1=tmask )
[2528]	281	!
[12377]	282	IF( ln_timing ) CALL timing_stop('ssh_atf')
[3294]	283	!
[12377]	284	END SUBROUTINE ssh_atf
[3]	285
[12377]	286	SUBROUTINE wAimp( kt, Kmm )
[10364]	287	!!----------------------------------------------------------------------
	288	!! * ROUTINE wAimp *
[12590]	289	!!
[10364]	290	!! ** Purpose : compute the Courant number and partition vertical velocity
	291	!! if a proportion needs to be treated implicitly
	292	!!
[12590]	293	!! ** Method : -
[10364]	294	!!
[12377]	295	!! ** action : ww : now vertical velocity (to be handled explicitly)
[10364]	296	!! : wi : now vertical velocity (for implicit treatment)
	297	!!
[11414]	298	!! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
[12590]	299	!! implicit scheme for vertical advection in oceanic modeling.
[11414]	300	!! Ocean Modelling, 91, 38-69.
[10364]	301	!!----------------------------------------------------------------------
	302	INTEGER, INTENT(in) :: kt ! time step
[12377]	303	INTEGER, INTENT(in) :: Kmm ! time level index
[10364]	304	!
	305	INTEGER :: ji, jj, jk ! dummy loop indices
[11293]	306	REAL(wp) :: zCu, zcff, z1_e3t ! local scalars
[10364]	307	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
[11407]	308	REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters
[10364]	309	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
	310	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
	311	!!----------------------------------------------------------------------
	312	!
	313	IF( ln_timing ) CALL timing_start('wAimp')
	314	!
	315	IF( kt == nit000 ) THEN
	316	IF(lwp) WRITE(numout,*)
	317	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
	318	IF(lwp) WRITE(numout,*) '~~~~~ '
[11293]	319	wi(:,:,:) = 0._wp
[10364]	320	ENDIF
	321	!
[11414]	322	! Calculate Courant numbers
	323	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
[12377]	324	DO_3D_00_00( 1, jpkm1 )
	325	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
	326	! 2*rdt and not r2dt (for restartability)
[12622]	327	Cu_adv(ji,jj,jk) = 2._wp * rdt * &
	328	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
	329	& + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) &
	330	& * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - &
	331	& MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) &
	332	& * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) &
	333	& * r1_e1e2t(ji ,jj) &
	334	& + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) &
	335	& * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - &
	336	& MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) &
	337	& * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) &
	338	& * r1_e1e2t(ji,jj ) &
	339	& ) * z1_e3t
[12377]	340	END_3D
[11414]	341	ELSE
[12377]	342	DO_3D_00_00( 1, jpkm1 )
	343	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
	344	! 2*rdt and not r2dt (for restartability)
[12622]	345	Cu_adv(ji,jj,jk) = 2._wp * rdt * &
	346	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
	347	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm), 0._wp ) - &
	348	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm), 0._wp ) ) &
	349	& * r1_e1e2t(ji,jj) &
	350	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm), 0._wp ) - &
	351	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm), 0._wp ) ) &
	352	& * r1_e1e2t(ji,jj) &
	353	& ) * z1_e3t
[12377]	354	END_3D
[11414]	355	ENDIF
[10907]	356	CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
[10364]	357	!
	358	CALL iom_put("Courant",Cu_adv)
	359	!
	360	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
[12377]	361	DO_3DS_11_11( jpkm1, 2, -1 )
	362	!
	363	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
[11293]	364	! alt:
[12590]	365	! IF ( ww(ji,jj,jk) > 0._wp ) THEN
	366	! zCu = Cu_adv(ji,jj,jk)
[11293]	367	! ELSE
	368	! zCu = Cu_adv(ji,jj,jk-1)
[12590]	369	! ENDIF
[12377]	370	!
	371	IF( zCu <= Cu_min ) THEN !<-- Fully explicit
	372	zcff = 0._wp
	373	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
	374	zcff = ( zCu - Cu_min )**2
	375	zcff = zcff / ( Fcu + zcff )
	376	ELSE !<-- Mostly implicit
	377	zcff = ( zCu - Cu_max )/ zCu
	378	ENDIF
	379	zcff = MIN(1._wp, zcff)
	380	!
	381	wi(ji,jj,jk) = zcff * ww(ji,jj,jk)
	382	ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
	383	!
	384	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl
	385	END_3D
[12590]	386	Cu_adv(:,:,1) = 0._wp
[10364]	387	ELSE
	388	! Fully explicit everywhere
[11407]	389	Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl
[10907]	390	wi (:,:,:) = 0._wp
[10364]	391	ENDIF
[12590]	392	CALL iom_put("wimp",wi)
[10364]	393	CALL iom_put("wi_cff",Cu_adv)
[12377]	394	CALL iom_put("wexp",ww)
[10364]	395	!
	396	IF( ln_timing ) CALL timing_stop('wAimp')
	397	!
	398	END SUBROUTINE wAimp
[3]	399	!!======================================================================
[1565]	400	END MODULE sshwzv

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

Download in other formats: