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.

traldf_bilapg.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 2787

Last change on this file since 2787 was 2715, checked in by rblod, 13 years ago
First attempt to put dynamic allocation on the trunk
Property svn:keywords set to `Id`
File size: 17.7 KB

Rev	Line
[3]	1	MODULE traldf_bilapg
	2	!!==============================================================================
	3	!! * MODULE traldf_bilapg *
[2528]	4	!! Ocean tracers: horizontal component of the lateral tracer mixing trend
[3]	5	!!==============================================================================
[2528]	6	!! History : 8.0 ! 1997-07 (G. Madec) Original code
	7	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	8	!! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC
	9	!!==============================================================================
[3]	10	#if defined key_ldfslp \|\| defined key_esopa
	11	!!----------------------------------------------------------------------
	12	!! 'key_ldfslp' rotation of the lateral mixing tensor
	13	!!----------------------------------------------------------------------
	14	!! tra_ldf_bilapg : update the tracer trend with the horizontal diffusion
	15	!! using an horizontal biharmonic operator in s-coordinate
	16	!! ldfght : ???
	17	!!----------------------------------------------------------------------
	18	USE oce ! ocean dynamics and tracers variables
	19	USE dom_oce ! ocean space and time domain variables
[74]	20	USE ldftra_oce ! ocean active tracers: lateral physics
[3]	21	USE in_out_manager ! I/O manager
	22	USE ldfslp ! iso-neutral slopes available
	23	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[132]	24	USE diaptr ! poleward transport diagnostics
[2528]	25	USE trc_oce ! share passive tracers/Ocean variables
[2715]	26	USE lib_mpp ! MPP library
[3]	27
	28	IMPLICIT NONE
	29	PRIVATE
	30
[2528]	31	PUBLIC tra_ldf_bilapg ! routine called by step.F90
[3]	32
	33	!! * Substitutions
	34	# include "domzgr_substitute.h90"
	35	# include "ldftra_substitute.h90"
	36	# include "ldfeiv_substitute.h90"
	37	!!----------------------------------------------------------------------
[2528]	38	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	39	!! $Id$
	40	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	41	!!----------------------------------------------------------------------
	42	CONTAINS
	43
[2528]	44	SUBROUTINE tra_ldf_bilapg( kt, cdtype, ptb, pta, kjpt )
[3]	45	!!----------------------------------------------------------------------
	46	!! * ROUTINE tra_ldf_bilapg *
	47	!!
[2528]	48	!! ** Purpose : Compute the before horizontal tracer diffusive
[3]	49	!! trend and add it to the general trend of tracer equation.
	50	!!
	51	!! ** Method : The lateral diffusive trends is provided by a 4th order
	52	!! operator rotated along geopotential surfaces. It is computed
	53	!! using before fields (forward in time) and geopotential slopes
	54	!! computed in routine inildf.
	55	!! -1- compute the geopotential harmonic operator applied to
[2528]	56	!! ptb and multiply it by the eddy diffusivity coefficient
	57	!! (done by a call to ldfght routine, result in wk1 arrays).
[3]	58	!! Applied the domain lateral boundary conditions by call to lbc_lnk
	59	!! -2- compute the geopotential harmonic operator applied to
[2528]	60	!! wk1 by a second call to ldfght routine (result in wk2)
[3]	61	!! arrays).
[2528]	62	!! -3- Add this trend to the general trend
	63	!! pta = pta + wk2
[3]	64	!!
[2528]	65	!! ** Action : - Update pta arrays with the before geopotential
[3]	66	!! biharmonic mixing trend.
[2528]	67	!!----------------------------------------------------------------------
[2715]	68	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	69	USE wrk_nemo, ONLY: wk1 => wrk_4d_1 , wk2 => wrk_4d_2 ! 4D workspace
	70	!
[2528]	71	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	72	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	73	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	74	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
	75	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2715]	76	!
	77	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[3]	78	!!----------------------------------------------------------------------
[132]	79
[2715]	80	IF( wrk_in_use(4, 1,2) ) THEN
	81	CALL ctl_stop('tra_ldf_bilapg: requested workspace arrays unavailable') ; RETURN
	82	ENDIF
	83
[2528]	84	IF( kt == nit000 ) THEN
[3]	85	IF(lwp) WRITE(numout,*)
[2528]	86	IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
[3]	87	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
	88	ENDIF
	89
[2528]	90	! 1. Laplacian of ptb * aht
[3]	91	! -----------------------------
[2528]	92	CALL ldfght( kt, cdtype, ptb, wk1, kjpt, 1 ) ! rotated harmonic operator applied to ptb and multiply by aht
	93	! ! output in wk1
	94	!
	95	DO jn = 1, kjpt
	96	CALL lbc_lnk( wk1(:,:,:,jn) , 'T', 1. ) ! Lateral boundary conditions on wk1 (unchanged sign)
	97	END DO
[3]	98
[2528]	99	! 2. Bilaplacian of ptb
[3]	100	! -------------------------
[2528]	101	CALL ldfght( kt, cdtype, wk1, wk2, kjpt, 2 ) ! rotated harmonic operator applied to wk1 ; output in wk2
[3]	102
	103
	104	! 3. Update the tracer trends (j-slab : 2, jpj-1)
	105	! ---------------------------
[2528]	106	DO jn = 1, kjpt
	107	DO jj = 2, jpjm1
	108	DO jk = 1, jpkm1
	109	DO ji = 2, jpim1
	110	! add it to the general tracer trends
	111	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + wk2(ji,jj,jk,jn)
	112	END DO
[3]	113	END DO
	114	END DO
[2528]	115	END DO
	116	!
[2715]	117	IF( wrk_not_released(4, 1,2) ) CALL ctl_stop('tra_ldf_bilapg : failed to release workspace arrays.')
	118	!
[3]	119	END SUBROUTINE tra_ldf_bilapg
	120
	121
[2528]	122	SUBROUTINE ldfght ( kt, cdtype, pt, plt, kjpt, kaht )
[3]	123	!!----------------------------------------------------------------------
	124	!! * ROUTINE ldfght *
	125	!!
[2528]	126	!! ** Purpose : Apply a geopotential harmonic operator to (pt) and
[3]	127	!! multiply it by the eddy diffusivity coefficient (if kaht=1).
	128	!! Routine only used in s-coordinates (l_sco=T) with bilaplacian
	129	!! operator (ln_traldf_bilap=T) acting along geopotential surfaces
	130	!! (ln_traldf_hor).
	131	!!
	132	!! ** Method : The harmonic operator rotated along geopotential
[2528]	133	!! surfaces is applied to (pt) using the slopes of geopotential
[3]	134	!! surfaces computed in inildf routine. The result is provided in
	135	!! (plt,pls) arrays. It is computed in 2 steps:
	136	!!
	137	!! First step: horizontal part of the operator. It is computed on
	138	!! ========== pt as follows (idem on ps)
	139	!! horizontal fluxes :
	140	!! zftu = e2ue3u/e1u di[ pt ] - e2uuslp dk[ mi(mk(pt)) ]
	141	!! zftv = e1ve3v/e2v dj[ pt ] - e1vvslp dk[ mj(mk(pt)) ]
	142	!! take the horizontal divergence of the fluxes (no divided by
	143	!! the volume element :
	144	!! plt = di-1[ zftu ] + dj-1[ zftv ]
	145	!!
	146	!! Second step: vertical part of the operator. It is computed on
	147	!! =========== pt as follows (idem on ps)
	148	!! vertical fluxes :
	149	!! zftw = e1te2t/e3w (wslpi^2+wslpj^2) dk-1[ pt ]
	150	!! - e2t * wslpi di[ mi(mk(pt)) ]
	151	!! - e1t * wslpj dj[ mj(mk(pt)) ]
	152	!! take the vertical divergence of the fluxes add it to the hori-
	153	!! zontal component, divide the result by the volume element and
	154	!! if kaht=1, multiply by the eddy diffusivity coefficient:
	155	!! plt = aht / (e1te2te3t) { plt + dk[ zftw ] }
	156	!! else:
	157	!! plt = 1 / (e1te2te3t) { plt + dk[ zftw ] }
	158	!!
[2528]	159	!!----------------------------------------------------------------------
[2715]	160	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, wrk_in_use_xz, wrk_not_released_xz
	161	USE oce , ONLY: zftv => ua ! ua used as workspace
	162	USE wrk_nemo, ONLY: zftu => wrk_2d_1 , zdkt => wrk_2d_2 , zdk1t => wrk_2d_3
	163	USE wrk_nemo, ONLY: zftw => wrk_xz_1 , zdit => wrk_xz_2
	164	USE wrk_nemo, ONLY: zdjt => wrk_xz_3 , zdj1t => wrk_xz_4
	165	!
[2528]	166	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	167	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	168	INTEGER , INTENT(in ) :: kjpt !: dimension of
	169	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: pt ! tracer fields ( before for 1st call
	170	! ! and laplacian of these fields for 2nd call.
	171	REAL(wp) , INTENT(out), DIMENSION(jpi,jpj,jpk,kjpt) :: plt !: partial harmonic operator applied to pt components except
	172	! !: second order vertical derivative term
	173	INTEGER , INTENT(in ) :: kaht !: =1 multiply the laplacian by the eddy diffusivity coeff.
	174	! !: =2 no multiplication
	175	!!
	176	INTEGER :: ji, jj, jk,jn ! dummy loop indices
	177	! ! temporary scalars
	178	REAL(wp) :: zabe1, zabe2, zmku, zmkv
	179	REAL(wp) :: zbtr, ztah, ztav
	180	REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4
[3]	181	!!----------------------------------------------------------------------
	182
[2715]	183	IF( wrk_in_use(2, 1,2,3) .OR. wrk_in_use_xz(1,2,3,4) )THEN
	184	CALL ctl_stop('ldfght : requested workspace arrays unavailable') ; RETURN
	185	ENDIF
[2528]	186	!
	187	DO jn = 1, kjpt
	188	! ! ********** ! ! ===============
	189	DO jk = 1, jpkm1 ! First step ! ! Horizontal slab
	190	! ! ********** ! ! ===============
	191
	192	! I.1 Vertical gradient of pt and ps at level jk and jk+1
	193	! -------------------------------------------------------
	194	! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
	195
	196	zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
	197	IF( jk == 1 ) THEN
	198	zdkt(:,:) = zdk1t(:,:)
	199	ELSE
	200	zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
	201	ENDIF
[3]	202
	203
[2528]	204	! I.2 Horizontal fluxes
	205	! ---------------------
	206
	207	DO jj = 1, jpjm1
	208	DO ji = 1, jpim1
	209	zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
	210	zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
	211
	212	zmku = 1./MAX( tmask(ji+1,jj,jk )+tmask(ji,jj,jk+1) &
	213	& +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk ),1. )
	214	zmkv = 1./MAX( tmask(ji,jj+1,jk )+tmask(ji,jj,jk+1) &
	215	& +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk ),1. )
	216
	217	zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
	218	zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
	219
	220	zftu(ji,jj)= umask(ji,jj,jk) * &
	221	& ( zabe1 *( pt (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) ) &
	222	& + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj) &
	223	& +zdk1t(ji+1,jj) + zdkt (ji,jj) ) )
	224
	225	zftv(ji,jj,jk)= vmask(ji,jj,jk) * &
	226	& ( zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) &
	227	& + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj) &
	228	& +zdk1t(ji,jj+1) + zdkt (ji,jj) ) )
	229	END DO
[3]	230	END DO
	231
	232
[2528]	233	! I.3 Second derivative (divergence) (not divided by the volume)
	234	! ---------------------
	235
	236	DO jj = 2 , jpjm1
	237	DO ji = 2 , jpim1
	238	ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
	239	plt(ji,jj,jk,jn) = ztah
	240	END DO
[3]	241	END DO
[2528]	242	! ! ===============
	243	END DO ! End of slab
	244	! ! ===============
	245	! "Poleward" diffusive heat or salt transport
	246	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
	247	IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) )
	248	IF( jn == jp_sal) str_ldf(:) = ptr_vj( zftv(:,:,:) )
	249	ENDIF
[3]	250
[2528]	251	! ! ************ ! ! ===============
	252	DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab
	253	! ! ************ ! ! ===============
	254
	255	! II.1 horizontal tracer gradient
	256	! -------------------------------
	257
	258	DO jk = 1, jpk
	259	DO ji = 1, jpim1
	260	zdit (ji,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj ,jk,jn) ) * umask(ji,jj ,jk)
	261	zdjt (ji,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj ,jk,jn) ) * vmask(ji,jj ,jk)
	262	zdj1t(ji,jk) = ( pt(ji ,jj ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
	263	END DO
[3]	264	END DO
	265
	266
[2528]	267	! II.2 Vertical fluxes
	268	! --------------------
	269
	270	! Surface and bottom vertical fluxes set to zero
	271	zftw(:, 1 ) = 0.e0
	272	zftw(:,jpk) = 0.e0
	273
	274	! interior (2=<jk=<jpk-1)
	275	DO jk = 2, jpkm1
	276	DO ji = 2, jpim1
	277	zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk) &
	278	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	279	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
	280
	281	zmku = 1./MAX( umask(ji ,jj,jk-1)+umask(ji-1,jj,jk) &
	282	& +umask(ji-1,jj,jk-1)+umask(ji ,jj,jk), 1. )
	283
	284	zmkv = 1./MAX( vmask(ji,jj ,jk-1)+vmask(ji,jj-1,jk) &
	285	& +vmask(ji,jj-1,jk-1)+vmask(ji,jj ,jk), 1. )
	286
	287	zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
	288	zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
	289
	290	zftw(ji,jk) = tmask(ji,jj,jk) * &
	291	& ( zcof0 * ( pt (ji,jj,jk-1,jn) - pt (ji ,jj,jk,jn) ) &
	292	& + zcof3 * ( zdit (ji ,jk-1 ) + zdit (ji-1,jk ) &
	293	& +zdit (ji-1 ,jk-1 ) + zdit (ji ,jk ) ) &
	294	& + zcof4 * ( zdjt (ji ,jk-1 ) + zdj1t(ji ,jk) &
	295	& +zdj1t(ji ,jk-1 ) + zdjt (ji ,jk ) ) )
	296	END DO
[3]	297	END DO
	298
	299
[2528]	300	! II.3 Divergence of vertical fluxes added to the horizontal divergence
	301	! ---------------------------------------------------------------------
	302
	303	IF( kaht == 1 ) THEN
	304	! multiply the laplacian by the eddy diffusivity coefficient
	305	DO jk = 1, jpkm1
	306	DO ji = 2, jpim1
	307	! eddy coef. divided by the volume element
	308	zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	309	! vertical divergence
	310	ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
	311	! harmonic operator applied to (pt,ps) and multiply by aht
	312	plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
	313	END DO
[3]	314	END DO
[2528]	315	ELSEIF( kaht == 2 ) THEN
	316	! second call, no multiplication
	317	DO jk = 1, jpkm1
	318	DO ji = 2, jpim1
	319	! inverse of the volume element
	320	zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	321	! vertical divergence
	322	ztav = zftw(ji,jk) - zftw(ji,jk+1)
	323	! harmonic operator applied to (pt,ps)
	324	plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
	325	END DO
[3]	326	END DO
[2528]	327	ELSE
	328	IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
	329	IF(lwp) WRITE(numout,*) ' We stop'
	330	STOP 'ldfght'
	331	ENDIF
	332	! ! ===============
	333	END DO ! End of slab
	334	! ! ===============
	335	END DO
	336	!
[2715]	337	IF( wrk_not_released(2, 1,2,3) .OR. &
	338	wrk_not_released_xz(1,2,3,4) ) CALL ctl_stop('ldfght : failed to release workspace arrays.')
	339	!
[3]	340	END SUBROUTINE ldfght
	341
	342	#else
	343	!!----------------------------------------------------------------------
	344	!! Dummy module : NO rotation of the lateral mixing tensor
	345	!!----------------------------------------------------------------------
	346	CONTAINS
[2528]	347	SUBROUTINE tra_ldf_bilapg( kt, cdtype, ptb, pta, kjpt ) ! Empty routine
	348	CHARACTER(len=3) :: cdtype
	349	REAL, DIMENSION(:,:,:,:) :: ptb, pta
	350	WRITE(,) 'tra_ldf_iso: You should not have seen this print! error?', kt, cdtype, ptb(1,1,1,1), pta(1,1,1,1), kjpt
[3]	351	END SUBROUTINE tra_ldf_bilapg
	352	#endif
	353
	354	!!==============================================================================
	355	END MODULE traldf_bilapg

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

Download in other formats: