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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 4401

Last change on this file since 4401 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

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

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

Download in other formats: