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.

dynzdf_imp.F90 in branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 5989

Last change on this file since 5989 was 5989, checked in by deazer, 8 years ago
Merging TMB and 25h diagnostics to head of trunk added brief documentation
Property svn:keywords set to `Id`
File size: 16.1 KB

Rev	Line
[3]	1	MODULE dynzdf_imp
[2715]	2	!!======================================================================
[3]	3	!! * MODULE dynzdf_imp *
	4	!! Ocean dynamics: vertical component(s) of the momentum mixing trend
[2715]	5	!!======================================================================
[2528]	6	!! History : OPA ! 1990-10 (B. Blanke) Original code
	7	!! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal
[2715]	8	!! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module
[2528]	9	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps
[3294]	10	!! 3.4 ! 2012-01 (H. Liu) Semi-implicit bottom friction
[503]	11	!!----------------------------------------------------------------------
[3]	12
	13	!!----------------------------------------------------------------------
[2715]	14	!! dyn_zdf_imp : update the momentum trend with the vertical diffusion using a implicit time-stepping
[3]	15	!!----------------------------------------------------------------------
	16	USE oce ! ocean dynamics and tracers
	17	USE dom_oce ! ocean space and time domain
[4292]	18	USE domvvl ! variable volume
[888]	19	USE sbc_oce ! surface boundary condition: ocean
	20	USE zdf_oce ! ocean vertical physics
[719]	21	USE phycst ! physical constants
[3]	22	USE in_out_manager ! I/O manager
[2715]	23	USE lib_mpp ! MPP library
[3294]	24	USE zdfbfr ! Bottom friction setup
	25	USE wrk_nemo ! Memory Allocation
	26	USE timing ! Timing
[5989]	27	USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form
[3]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[2528]	32	PUBLIC dyn_zdf_imp ! called by step.F90
[3]	33
[4292]	34	REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
	35
[3]	36	!! * Substitutions
	37	# include "domzgr_substitute.h90"
	38	# include "vectopt_loop_substitute.h90"
	39	!!----------------------------------------------------------------------
[2528]	40	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	41	!! $Id$
[2528]	42	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	43	!!----------------------------------------------------------------------
	44	CONTAINS
	45
[503]	46	SUBROUTINE dyn_zdf_imp( kt, p2dt )
[3]	47	!!----------------------------------------------------------------------
	48	!! * ROUTINE dyn_zdf_imp *
	49	!!
	50	!! ** Purpose : Compute the trend due to the vert. momentum diffusion
	51	!! and the surface forcing, and add it to the general trend of
	52	!! the momentum equations.
	53	!!
	54	!! ** Method : The vertical momentum mixing trend is given by :
	55	!! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
	56	!! backward time stepping
[2528]	57	!! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
[3]	58	!! Bottom boundary conditions : bottom stress (cf zdfbfr.F)
	59	!! Add this trend to the general trend ua :
	60	!! ua = ua + dz( avmu dz(u) )
	61	!!
[2528]	62	!! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
[3]	63	!!---------------------------------------------------------------------
[3294]	64	INTEGER , INTENT(in) :: kt ! ocean time-step index
[2715]	65	REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step
[2528]	66	!!
[2715]	67	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	68	INTEGER :: ikbu, ikbv ! local integers
[2715]	69	REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars
[4292]	70	REAL(wp) :: ze3ua, ze3va
[3294]	71	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws
	72	!!----------------------------------------------------------------------
	73	!
	74	IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_imp')
	75	!
	76	CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws )
	77	!
[3]	78	IF( kt == nit000 ) THEN
	79	IF(lwp) WRITE(numout,*)
	80	IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
	81	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
[4292]	82	!
	83	IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator
	84	ELSE ; r_vvl = 0._wp
	85	ENDIF
[3]	86	ENDIF
	87
[5989]	88	! 1. Time step dynamics
	89	! ---------------------
[455]	90
[5989]	91	IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity
	92	DO jk = 1, jpkm1
	93	ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
	94	va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
	95	END DO
	96	ELSE ! applied on thickness weighted velocity
	97	DO jk = 1, jpkm1
	98	ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &
	99	& + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) &
	100	& / fse3u_a(:,:,jk) * umask(:,:,jk)
	101	va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &
	102	& + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) &
	103	& / fse3v_a(:,:,jk) * vmask(:,:,jk)
	104	END DO
	105	ENDIF
	106
	107	! 2. Apply semi-implicit bottom friction
[3294]	108	! --------------------------------------
	109	! Only needed for semi-implicit bottom friction setup. The explicit
	110	! bottom friction has been included in "u(v)a" which act as the R.H.S
	111	! column vector of the tri-diagonal matrix equation
	112	!
	113	IF( ln_bfrimp ) THEN
[4292]	114	DO jj = 2, jpjm1
	115	DO ji = 2, jpim1
	116	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
	117	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
	118	avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
	119	avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
	120	END DO
[3294]	121	END DO
[5260]	122	IF ( ln_isfcav ) THEN
	123	DO jj = 2, jpjm1
	124	DO ji = 2, jpim1
	125	ikbu = miku(ji,jj) ! ocean top level at u- and v-points
	126	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
	127	IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
	128	IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
	129	END DO
	130	END DO
	131	END IF
[3294]	132	ENDIF
	133
[5989]	134	! With split-explicit free surface, barotropic stress is treated explicitly
	135	! Update velocities at the bottom.
	136	! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
	137	! not lead to the effective stress seen over the whole barotropic loop.
	138	IF ( ln_bfrimp .AND.ln_dynspg_ts ) THEN
[4292]	139	! remove barotropic velocities:
	140	DO jk = 1, jpkm1
	141	ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
	142	va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
[5260]	143	END DO
	144	! Add bottom/top stress due to barotropic component only:
[4292]	145	DO jj = 2, jpjm1
	146	DO ji = fs_2, fs_jpim1 ! vector opt.
	147	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
	148	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
	149	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)
	150	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)
	151	ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
	152	va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
	153	END DO
	154	END DO
[5260]	155	IF ( ln_isfcav ) THEN
	156	DO jj = 2, jpjm1
	157	DO ji = fs_2, fs_jpim1 ! vector opt.
	158	ikbu = miku(ji,jj) ! top ocean level at u- and v-points
	159	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
	160	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)
	161	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)
	162	ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
	163	va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
	164	END DO
	165	END DO
	166	END IF
[4292]	167	ENDIF
	168
[5989]	169	! 3. Vertical diffusion on u
[3]	170	! ---------------------------
	171	! Matrix and second member construction
[1662]	172	! bottom boundary condition: both zwi and zws must be masked as avmu can take
[3294]	173	! non zero value at the ocean bottom depending on the bottom friction used.
[2528]	174	!
	175	DO jk = 1, jpkm1 ! Matrix
[3]	176	DO jj = 2, jpjm1
	177	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	178	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl * fse3u_a(ji,jj,jk) ! after scale factor at T-point
	179	zcoef = - p2dt / ze3ua
[5260]	180	zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )
	181	zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk )
	182	zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
	183	zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1)
	184	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
[3]	185	END DO
	186	END DO
	187	END DO
[4292]	188	DO jj = 2, jpjm1 ! Surface boundary conditions
[3]	189	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	190	zwi(ji,jj,1) = 0._wp
	191	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]	192	END DO
	193	END DO
	194
	195	! Matrix inversion starting from the first level
	196	!-----------------------------------------------------------------------
	197	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	198	!
	199	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	200	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	201	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	202	! ( ... )( ... ) ( ... )
	203	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	204	!
	205	! m is decomposed in the product of an upper and a lower triangular matrix
	206	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	207	! The solution (the after velocity) is in ua
	208	!-----------------------------------------------------------------------
[2528]	209	!
[5989]	210	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	211	DO jj = 2, jpjm1
	212	DO ji = fs_2, fs_jpim1 ! vector opt.
	213	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	214	END DO
	215	END DO
	216	END DO
[2528]	217	!
	218	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	219	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	220	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1)
	221	ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
[5260]	222	& / ( ze3ua * rau0 ) * umask(ji,jj,1)
[3]	223	END DO
	224	END DO
	225	DO jk = 2, jpkm1
[5260]	226	DO jj = 2, jpjm1
	227	DO ji = fs_2, fs_jpim1
[4292]	228	zrhs = ua(ji,jj,jk) ! zrhs=right hand side
[3]	229	ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
	230	END DO
	231	END DO
	232	END DO
[2528]	233	!
	234	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==
[3]	235	DO ji = fs_2, fs_jpim1 ! vector opt.
	236	ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	237	END DO
	238	END DO
	239	DO jk = jpk-2, 1, -1
[5260]	240	DO jj = 2, jpjm1
	241	DO ji = fs_2, fs_jpim1
[2528]	242	ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]	243	END DO
	244	END DO
	245	END DO
	246
[5989]	247	! 4. Vertical diffusion on v
[3]	248	! ---------------------------
	249	! Matrix and second member construction
[1662]	250	! bottom boundary condition: both zwi and zws must be masked as avmv can take
[3294]	251	! non zero value at the ocean bottom depending on the bottom friction used
[2528]	252	!
	253	DO jk = 1, jpkm1 ! Matrix
[3]	254	DO jj = 2, jpjm1
	255	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	256	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk) + r_vvl * fse3v_a(ji,jj,jk) ! after scale factor at T-point
	257	zcoef = - p2dt / ze3va
[2528]	258	zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )
[5260]	259	zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk)
[2528]	260	zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
[5260]	261	zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
	262	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
[3]	263	END DO
	264	END DO
	265	END DO
[4292]	266	DO jj = 2, jpjm1 ! Surface boundary conditions
[3]	267	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	268	zwi(ji,jj,1) = 0._wp
	269	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]	270	END DO
	271	END DO
	272
	273	! Matrix inversion
	274	!-----------------------------------------------------------------------
	275	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	276	!
	277	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	278	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	279	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	280	! ( ... )( ... ) ( ... )
	281	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	282	!
[2528]	283	! m is decomposed in the product of an upper and lower triangular matrix
[3]	284	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	285	! The solution (after velocity) is in 2d array va
	286	!-----------------------------------------------------------------------
[2528]	287	!
[5989]	288	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	289	DO jj = 2, jpjm1
	290	DO ji = fs_2, fs_jpim1 ! vector opt.
	291	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	292	END DO
	293	END DO
	294	END DO
[2528]	295	!
	296	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[5989]	297	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	298	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1)
	299	va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
	300	& / ( ze3va * rau0 )
[3]	301	END DO
	302	END DO
	303	DO jk = 2, jpkm1
	304	DO jj = 2, jpjm1
	305	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	306	zrhs = va(ji,jj,jk) ! zrhs=right hand side
[3]	307	va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
	308	END DO
	309	END DO
	310	END DO
[2528]	311	!
[4292]	312	DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==
[3]	313	DO ji = fs_2, fs_jpim1 ! vector opt.
	314	va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	315	END DO
	316	END DO
	317	DO jk = jpk-2, 1, -1
[5260]	318	DO jj = 2, jpjm1
	319	DO ji = fs_2, fs_jpim1
[2528]	320	va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]	321	END DO
	322	END DO
	323	END DO
	324
[4292]	325	! J. Chanut: Lines below are useless ?
[5989]	326	!! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces)
[3294]	327	IF( ln_bfrimp ) THEN
[5260]	328	DO jj = 2, jpjm1
	329	DO ji = 2, jpim1
	330	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
	331	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
	332	avmu(ji,jj,ikbu+1) = 0.e0
	333	avmv(ji,jj,ikbv+1) = 0.e0
	334	END DO
	335	END DO
	336	IF (ln_isfcav) THEN
	337	DO jj = 2, jpjm1
	338	DO ji = 2, jpim1
	339	ikbu = miku(ji,jj) ! ocean top level at u- and v-points
	340	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
	341	IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
	342	IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
	343	END DO
	344	END DO
	345	END IF
[3294]	346	ENDIF
[2528]	347	!
[3294]	348	CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)
[2715]	349	!
[3294]	350	IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp')
	351	!
[3]	352	END SUBROUTINE dyn_zdf_imp
	353
	354	!!==============================================================================
	355	END MODULE dynzdf_imp

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

Download in other formats: