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.

dynldf_bilap.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90 @ 2590

Last change on this file since 2590 was 2590, checked in by trackstand2, 13 years ago
Merge branch 'dynamic_memory' into master-svn-dyn
Property svn:keywords set to `Id`
File size: 9.7 KB

Rev	Line
[3]	1	MODULE dynldf_bilap
	2	!!======================================================================
	3	!! * MODULE dynldf_bilap *
	4	!! Ocean dynamics: lateral viscosity trend
	5	!!======================================================================
	6
	7	!!----------------------------------------------------------------------
	8	!! dyn_ldf_bilap : update the momentum trend with the lateral diffusion
	9	!! using an iso-level bilaplacian operator
	10	!!----------------------------------------------------------------------
	11	!! * Modules used
	12	USE oce ! ocean dynamics and tracers
	13	USE dom_oce ! ocean space and time domain
	14	USE ldfdyn_oce ! ocean dynamics: lateral physics
	15	USE in_out_manager ! I/O manager
[216]	16	USE trdmod ! ocean dynamics trends
	17	USE trdmod_oce ! ocean variables trends
[3]	18	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	19
	20	IMPLICIT NONE
	21	PRIVATE
	22
	23	!! * Routine accessibility
	24	PUBLIC dyn_ldf_bilap ! called by step.F90
	25
	26	!! * Substitutions
	27	# include "domzgr_substitute.h90"
	28	# include "ldfdyn_substitute.h90"
	29	# include "vectopt_loop_substitute.h90"
	30	!!----------------------------------------------------------------------
[2528]	31	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	32	!! $Id$
[2528]	33	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	34	!!----------------------------------------------------------------------
	35
	36	CONTAINS
	37
	38	SUBROUTINE dyn_ldf_bilap( kt )
	39	!!----------------------------------------------------------------------
	40	!! * ROUTINE dyn_ldf_bilap *
	41	!!
	42	!! ** Purpose : Compute the before trend of the lateral momentum
	43	!! diffusion and add it to the general trend of momentum equation.
	44	!!
	45	!! ** Method : The before horizontal momentum diffusion trend is a
	46	!! bi-harmonic operator (bilaplacian type) which separates the
	47	!! divergent and rotational parts of the flow.
	48	!! Its horizontal components are computed as follow:
	49	!! laplacian:
	50	!! zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ]
	51	!! zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ]
	52	!! third derivative:
	53	!! * multiply by the eddy viscosity coef. at u-, v-point, resp.
	54	!! zlu = ahmu * zlu
	55	!! zlv = ahmv * zlv
	56	!! * curl and divergence of the laplacian
	57	!! zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] )
	58	!! zut = 1/(e1te2te3t) ( di[e2ue3u zlu] + dj[e1ve3v zlv] )
	59	!! bilaplacian:
	60	!! diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ]
	61	!! diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ]
[455]	62	!! If ln_sco=F and ln_zps=F, the vertical scale factors in the
[3]	63	!! rotational part of the diffusion are simplified
	64	!! Add this before trend to the general trend (ua,va):
	65	!! (ua,va) = (ua,va) + (diffu,diffv)
	66	!! 'key_trddyn' defined: the two components of the horizontal
	67	!! diffusion trend are saved.
	68	!!
	69	!! ** Action : - Update (ua,va) with the before iso-level biharmonic
	70	!! mixing trend.
	71	!!
	72	!! History :
	73	!! ! 90-09 (G. Madec) Original code
	74	!! ! 91-11 (G. Madec)
	75	!! ! 93-03 (M. Guyon) symetrical conditions (M. Guyon)
	76	!! ! 96-01 (G. Madec) statement function for e3
	77	!! ! 97-07 (G. Madec) lbc calls
	78	!! 8.5 ! 02-08 (G. Madec) F90: Free form and module
[216]	79	!! 9.0 ! 04-08 (C. Talandier) New trends organization
[3]	80	!!----------------------------------------------------------------------
[2590]	81	USE wrk_nemo, ONLY: wrk_use, wrk_release
	82	USE wrk_nemo, ONLY: zcu => wrk_2d_1, zcv => wrk_2d_2
	83	USE wrk_nemo, ONLY: zuf => wrk_3d_1, zut => wrk_3d_2, &
	84	zlu => wrk_3d_3, zlv => wrk_3d_4
[3]	85	!! * Arguments
	86	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	87
	88	!! * Local declarations
	89	INTEGER :: ji, jj, jk ! dummy loop indices
	90	REAL(wp) :: zua, zva, zbt, ze2u, ze2v ! temporary scalar
	91	!!----------------------------------------------------------------------
	92	!! OPA 8.5, LODYC-IPSL (2002)
	93	!!----------------------------------------------------------------------
	94
[2590]	95	IF( (.NOT. wrk_use(2, 1,2)) .OR. (.NOT. wrk_use(3, 1,2,3,4)) )THEN
	96	CALL ctl_stop('dyn_ldf_bilap : requested workspace arrays unavailable.')
	97	RETURN
	98	END IF
	99
[3]	100	IF( kt == nit000 ) THEN
	101	IF(lwp) WRITE(numout,*)
	102	IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator'
	103	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
	104	ENDIF
[216]	105
[474]	106	!!bug gm this should be enough
	107	!!$ zuf(:,:,jpk) = 0.e0
	108	!!$ zut(:,:,jpk) = 0.e0
	109	!!$ zlu(:,:,jpk) = 0.e0
	110	!!$ zlv(:,:,jpk) = 0.e0
	111	zuf(:,:,:) = 0.e0
	112	zut(:,:,:) = 0.e0
	113	zlu(:,:,:) = 0.e0
	114	zlv(:,:,:) = 0.e0
	115
[3]	116	! ! ===============
	117	DO jk = 1, jpkm1 ! Horizontal slab
	118	! ! ===============
	119	! Laplacian
	120	! ---------
	121
[455]	122	IF( ln_sco .OR. ln_zps ) THEN ! s-coordinate or z-coordinate with partial steps
[474]	123	zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
[3]	124	DO jj = 2, jpjm1
	125	DO ji = fs_2, fs_jpim1 ! vector opt.
[474]	126	zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) &
[3]	127	& + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
	128
[474]	129	zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) &
[3]	130	& + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
	131	END DO
	132	END DO
[455]	133	ELSE ! z-coordinate - full step
[3]	134	DO jj = 2, jpjm1
	135	DO ji = fs_2, fs_jpim1 ! vector opt.
[474]	136	zlu(ji,jj,jk) = - ( rotb (ji ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj) &
[3]	137	& + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj ,jk) ) / e1u(ji,jj)
	138
[474]	139	zlv(ji,jj,jk) = + ( rotb (ji,jj ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj) &
[3]	140	& + ( hdivb(ji,jj+1,jk) - hdivb(ji ,jj,jk) ) / e2v(ji,jj)
	141	END DO
	142	END DO
	143	ENDIF
[474]	144	ENDDO
[3]	145
[474]	146	! Boundary conditions on the laplacian (zlu,zlv)
	147	CALL lbc_lnk( zlu, 'U', -1. )
	148	CALL lbc_lnk( zlv, 'V', -1. )
[3]	149
	150
[474]	151	DO jk = 1, jpkm1
	152
[3]	153	! Third derivative
	154	! ----------------
	155
	156	! Multiply by the eddy viscosity coef. (at u- and v-points)
[474]	157	zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk)
	158	zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk)
[3]	159
	160	! Contravariant "laplacian"
[474]	161	zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
	162	zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
[3]	163
	164	! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
	165	DO jj = 1, jpjm1
	166	DO ji = 1, fs_jpim1 ! vector opt.
[474]	167	zuf(ji,jj,jk) = fmask(ji,jj,jk) * ( zcv(ji+1,jj ) - zcv(ji,jj) &
[3]	168	& - zcu(ji ,jj+1) + zcu(ji,jj) ) &
	169	& * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
	170	END DO
	171	END DO
	172
	173	! Laplacian Horizontal fluxes
	174	DO jj = 1, jpjm1
	175	DO ji = 1, fs_jpim1 ! vector opt.
[474]	176	zlu(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj,jk)
	177	zlv(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj,jk)
[3]	178	END DO
	179	END DO
	180
	181	! Laplacian divergence
	182	DO jj = 2, jpj
	183	DO ji = fs_2, jpi ! vector opt.
	184	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
[474]	185	zut(ji,jj,jk) = ( zlu(ji,jj,jk) - zlu(ji-1,jj ,jk) &
	186	& + zlv(ji,jj,jk) - zlv(ji ,jj-1,jk) ) / zbt
[3]	187	END DO
	188	END DO
[474]	189	END DO
[3]	190
	191
[235]	192	! boundary conditions on the laplacian curl and div (zuf,zut)
[474]	193	!!bug gm no need to do this 2 following lbc...
[235]	194	CALL lbc_lnk( zuf, 'F', 1. )
	195	CALL lbc_lnk( zut, 'T', 1. )
[3]	196
[474]	197	DO jk = 1, jpkm1
	198
[3]	199	! Bilaplacian
	200	! -----------
	201
	202	DO jj = 2, jpjm1
	203	DO ji = fs_2, fs_jpim1 ! vector opt.
	204	ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
	205	ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
	206	! horizontal biharmonic diffusive trends
[474]	207	zua = - ( zuf(ji ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u &
	208	& + ( zut(ji+1,jj,jk) - zut(ji,jj ,jk) ) / e1u(ji,jj)
[3]	209
[474]	210	zva = + ( zuf(ji,jj ,jk) - zuf(ji-1,jj,jk) ) / ze2v &
	211	& + ( zut(ji,jj+1,jk) - zut(ji ,jj,jk) ) / e2v(ji,jj)
[3]	212	! add it to the general momentum trends
	213	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
	214	va(ji,jj,jk) = va(ji,jj,jk) + zva
	215	END DO
	216	END DO
	217
	218	! ! ===============
	219	END DO ! End of slab
	220	! ! ===============
[2590]	221	IF( (.NOT. wrk_release(2, 1,2)) .OR. &
	222	(.NOT. wrk_release(3, 1,2,3,4)) )THEN
	223	CALL ctl_stop('dyn_ldf_bilap : failed to release workspace arrays.')
	224	END IF
	225	!
[3]	226	END SUBROUTINE dyn_ldf_bilap
	227
	228	!!======================================================================
	229	END MODULE dynldf_bilap

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

Download in other formats: