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/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 5279

Last change on this file since 5279 was 5279, checked in by dancopsey, 9 years ago
Removed SVN keywords.
File size: 17.8 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
[3294]	27	USE wrk_nemo ! Memory Allocation
	28	USE timing ! Timing
[3]	29
	30	IMPLICIT NONE
	31	PRIVATE
	32
[2528]	33	PUBLIC tra_ldf_bilapg ! routine called by step.F90
[3]	34
	35	!! * Substitutions
	36	# include "domzgr_substitute.h90"
	37	# include "ldftra_substitute.h90"
	38	# include "ldfeiv_substitute.h90"
	39	!!----------------------------------------------------------------------
[2528]	40	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	41	!! $Id$
	42	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	43	!!----------------------------------------------------------------------
	44	CONTAINS
	45
[3294]	46	SUBROUTINE tra_ldf_bilapg( kt, kit000, cdtype, ptb, pta, kjpt )
[3]	47	!!----------------------------------------------------------------------
	48	!! * ROUTINE tra_ldf_bilapg *
	49	!!
[2528]	50	!! ** Purpose : Compute the before horizontal tracer diffusive
[3]	51	!! trend and add it to the general trend of tracer equation.
	52	!!
	53	!! ** Method : The lateral diffusive trends is provided by a 4th order
	54	!! operator rotated along geopotential surfaces. It is computed
	55	!! using before fields (forward in time) and geopotential slopes
	56	!! computed in routine inildf.
	57	!! -1- compute the geopotential harmonic operator applied to
[2528]	58	!! ptb and multiply it by the eddy diffusivity coefficient
	59	!! (done by a call to ldfght routine, result in wk1 arrays).
[3]	60	!! Applied the domain lateral boundary conditions by call to lbc_lnk
	61	!! -2- compute the geopotential harmonic operator applied to
[2528]	62	!! wk1 by a second call to ldfght routine (result in wk2)
[3]	63	!! arrays).
[2528]	64	!! -3- Add this trend to the general trend
	65	!! pta = pta + wk2
[3]	66	!!
[2528]	67	!! ** Action : - Update pta arrays with the before geopotential
[3]	68	!! biharmonic mixing trend.
[2528]	69	!!----------------------------------------------------------------------
[2715]	70	!
[2528]	71	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	72	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	73	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	74	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	75	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
	76	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2715]	77	!
	78	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[3294]	79	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zwk1, zwk2
[3]	80	!!----------------------------------------------------------------------
[3294]	81	!
	82	IF( nn_timing == 1 ) CALL timing_start('tra_ldf_bilapg')
	83	!
	84	CALL wrk_alloc( jpi, jpj, jpk, kjpt, zwk1, zwk2 )
	85	!
	86	IF( kt == kit000 ) THEN
[3]	87	IF(lwp) WRITE(numout,*)
[2528]	88	IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
[3]	89	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
	90	ENDIF
	91
[2528]	92	! 1. Laplacian of ptb * aht
[3]	93	! -----------------------------
[3294]	94	CALL ldfght( kt, cdtype, ptb, zwk1, kjpt, 1 ) ! rotated harmonic operator applied to ptb and multiply by aht
[2528]	95	! ! output in wk1
	96	!
	97	DO jn = 1, kjpt
[3294]	98	CALL lbc_lnk( zwk1(:,:,:,jn) , 'T', 1. ) ! Lateral boundary conditions on wk1 (unchanged sign)
[2528]	99	END DO
[3]	100
[2528]	101	! 2. Bilaplacian of ptb
[3]	102	! -------------------------
[3294]	103	CALL ldfght( kt, cdtype, zwk1, zwk2, kjpt, 2 ) ! rotated harmonic operator applied to wk1 ; output in wk2
[3]	104
	105
	106	! 3. Update the tracer trends (j-slab : 2, jpj-1)
	107	! ---------------------------
[2528]	108	DO jn = 1, kjpt
	109	DO jj = 2, jpjm1
	110	DO jk = 1, jpkm1
	111	DO ji = 2, jpim1
	112	! add it to the general tracer trends
[3294]	113	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zwk2(ji,jj,jk,jn)
[2528]	114	END DO
[3]	115	END DO
	116	END DO
[2528]	117	END DO
	118	!
[3294]	119	CALL wrk_dealloc( jpi, jpj, jpk, kjpt, zwk1, zwk2 )
[2715]	120	!
[3294]	121	IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_bilapg')
	122	!
[3]	123	END SUBROUTINE tra_ldf_bilapg
	124
	125
[2528]	126	SUBROUTINE ldfght ( kt, cdtype, pt, plt, kjpt, kaht )
[3]	127	!!----------------------------------------------------------------------
	128	!! * ROUTINE ldfght *
	129	!!
[2528]	130	!! ** Purpose : Apply a geopotential harmonic operator to (pt) and
[3]	131	!! multiply it by the eddy diffusivity coefficient (if kaht=1).
	132	!! Routine only used in s-coordinates (l_sco=T) with bilaplacian
	133	!! operator (ln_traldf_bilap=T) acting along geopotential surfaces
	134	!! (ln_traldf_hor).
	135	!!
	136	!! ** Method : The harmonic operator rotated along geopotential
[2528]	137	!! surfaces is applied to (pt) using the slopes of geopotential
[3]	138	!! surfaces computed in inildf routine. The result is provided in
	139	!! (plt,pls) arrays. It is computed in 2 steps:
	140	!!
	141	!! First step: horizontal part of the operator. It is computed on
	142	!! ========== pt as follows (idem on ps)
	143	!! horizontal fluxes :
	144	!! zftu = e2ue3u/e1u di[ pt ] - e2uuslp dk[ mi(mk(pt)) ]
	145	!! zftv = e1ve3v/e2v dj[ pt ] - e1vvslp dk[ mj(mk(pt)) ]
	146	!! take the horizontal divergence of the fluxes (no divided by
	147	!! the volume element :
	148	!! plt = di-1[ zftu ] + dj-1[ zftv ]
	149	!!
	150	!! Second step: vertical part of the operator. It is computed on
	151	!! =========== pt as follows (idem on ps)
	152	!! vertical fluxes :
	153	!! zftw = e1te2t/e3w (wslpi^2+wslpj^2) dk-1[ pt ]
	154	!! - e2t * wslpi di[ mi(mk(pt)) ]
	155	!! - e1t * wslpj dj[ mj(mk(pt)) ]
	156	!! take the vertical divergence of the fluxes add it to the hori-
	157	!! zontal component, divide the result by the volume element and
	158	!! if kaht=1, multiply by the eddy diffusivity coefficient:
	159	!! plt = aht / (e1te2te3t) { plt + dk[ zftw ] }
	160	!! else:
	161	!! plt = 1 / (e1te2te3t) { plt + dk[ zftw ] }
	162	!!
[2528]	163	!!----------------------------------------------------------------------
[2715]	164	USE oce , ONLY: zftv => ua ! ua used as workspace
	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
[3294]	181	REAL(wp), POINTER, DIMENSION(:,:) :: zftu, zdkt, zdk1t
	182	REAL(wp), POINTER, DIMENSION(:,:) :: zftw, zdit, zdjt, zdj1t
[3]	183	!!----------------------------------------------------------------------
[2528]	184	!
[3294]	185	IF( nn_timing == 1 ) CALL timing_start('ldfght')
	186	!
	187	CALL wrk_alloc( jpi, jpj, zftu, zdkt, zdk1t )
	188	CALL wrk_alloc( jpi, jpk, zftw, zdit, zdjt, zdj1t )
	189	!
[2528]	190	DO jn = 1, kjpt
	191	! ! ********** ! ! ===============
	192	DO jk = 1, jpkm1 ! First step ! ! Horizontal slab
	193	! ! ********** ! ! ===============
	194
	195	! I.1 Vertical gradient of pt and ps at level jk and jk+1
	196	! -------------------------------------------------------
	197	! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
	198
	199	zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
	200	IF( jk == 1 ) THEN
	201	zdkt(:,:) = zdk1t(:,:)
	202	ELSE
	203	zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
	204	ENDIF
[3]	205
	206
[2528]	207	! I.2 Horizontal fluxes
	208	! ---------------------
	209
	210	DO jj = 1, jpjm1
	211	DO ji = 1, jpim1
[4292]	212	zabe1 = re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
	213	zabe2 = re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
[2528]	214
	215	zmku = 1./MAX( tmask(ji+1,jj,jk )+tmask(ji,jj,jk+1) &
	216	& +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk ),1. )
	217	zmkv = 1./MAX( tmask(ji,jj+1,jk )+tmask(ji,jj,jk+1) &
	218	& +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk ),1. )
	219
	220	zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
	221	zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
	222
	223	zftu(ji,jj)= umask(ji,jj,jk) * &
	224	& ( zabe1 *( pt (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) ) &
	225	& + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj) &
	226	& +zdk1t(ji+1,jj) + zdkt (ji,jj) ) )
	227
	228	zftv(ji,jj,jk)= vmask(ji,jj,jk) * &
	229	& ( zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) &
	230	& + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj) &
	231	& +zdk1t(ji,jj+1) + zdkt (ji,jj) ) )
	232	END DO
[3]	233	END DO
	234
	235
[2528]	236	! I.3 Second derivative (divergence) (not divided by the volume)
	237	! ---------------------
	238
	239	DO jj = 2 , jpjm1
	240	DO ji = 2 , jpim1
	241	ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
	242	plt(ji,jj,jk,jn) = ztah
	243	END DO
[3]	244	END DO
[2528]	245	! ! ===============
	246	END DO ! End of slab
	247	! ! ===============
	248	! "Poleward" diffusive heat or salt transport
	249	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
[3805]	250	! note sign is reversed to give down-gradient diffusive transports (#1043)
	251	IF( jn == jp_tem) htr_ldf(:) = ptr_vj( -zftv(:,:,:) )
	252	IF( jn == jp_sal) str_ldf(:) = ptr_vj( -zftv(:,:,:) )
[2528]	253	ENDIF
[3]	254
[2528]	255	! ! ************ ! ! ===============
	256	DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab
	257	! ! ************ ! ! ===============
	258
	259	! II.1 horizontal tracer gradient
	260	! -------------------------------
	261
	262	DO jk = 1, jpk
	263	DO ji = 1, jpim1
	264	zdit (ji,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj ,jk,jn) ) * umask(ji,jj ,jk)
	265	zdjt (ji,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj ,jk,jn) ) * vmask(ji,jj ,jk)
	266	zdj1t(ji,jk) = ( pt(ji ,jj ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
	267	END DO
[3]	268	END DO
	269
	270
[2528]	271	! II.2 Vertical fluxes
	272	! --------------------
	273
	274	! Surface and bottom vertical fluxes set to zero
	275	zftw(:, 1 ) = 0.e0
	276	zftw(:,jpk) = 0.e0
	277
	278	! interior (2=<jk=<jpk-1)
	279	DO jk = 2, jpkm1
	280	DO ji = 2, jpim1
[4292]	281	zcof0 = e12t(ji,jj) / fse3w_n(ji,jj,jk) &
[2528]	282	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	283	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
	284
	285	zmku = 1./MAX( umask(ji ,jj,jk-1)+umask(ji-1,jj,jk) &
	286	& +umask(ji-1,jj,jk-1)+umask(ji ,jj,jk), 1. )
	287
	288	zmkv = 1./MAX( vmask(ji,jj ,jk-1)+vmask(ji,jj-1,jk) &
	289	& +vmask(ji,jj-1,jk-1)+vmask(ji,jj ,jk), 1. )
	290
	291	zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
	292	zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
	293
	294	zftw(ji,jk) = tmask(ji,jj,jk) * &
	295	& ( zcof0 * ( pt (ji,jj,jk-1,jn) - pt (ji ,jj,jk,jn) ) &
	296	& + zcof3 * ( zdit (ji ,jk-1 ) + zdit (ji-1,jk ) &
	297	& +zdit (ji-1 ,jk-1 ) + zdit (ji ,jk ) ) &
	298	& + zcof4 * ( zdjt (ji ,jk-1 ) + zdj1t(ji ,jk) &
	299	& +zdj1t(ji ,jk-1 ) + zdjt (ji ,jk ) ) )
	300	END DO
[3]	301	END DO
	302
	303
[2528]	304	! II.3 Divergence of vertical fluxes added to the horizontal divergence
	305	! ---------------------------------------------------------------------
	306
	307	IF( kaht == 1 ) THEN
	308	! multiply the laplacian by the eddy diffusivity coefficient
	309	DO jk = 1, jpkm1
	310	DO ji = 2, jpim1
	311	! eddy coef. divided by the volume element
[4292]	312	zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
[2528]	313	! vertical divergence
	314	ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
	315	! harmonic operator applied to (pt,ps) and multiply by aht
	316	plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
	317	END DO
[3]	318	END DO
[2528]	319	ELSEIF( kaht == 2 ) THEN
	320	! second call, no multiplication
	321	DO jk = 1, jpkm1
	322	DO ji = 2, jpim1
	323	! inverse of the volume element
[4292]	324	zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
[2528]	325	! vertical divergence
	326	ztav = zftw(ji,jk) - zftw(ji,jk+1)
	327	! harmonic operator applied to (pt,ps)
	328	plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
	329	END DO
[3]	330	END DO
[2528]	331	ELSE
	332	IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
	333	IF(lwp) WRITE(numout,*) ' We stop'
	334	STOP 'ldfght'
	335	ENDIF
	336	! ! ===============
	337	END DO ! End of slab
	338	! ! ===============
	339	END DO
	340	!
[3294]	341	CALL wrk_dealloc( jpi, jpj, zftu, zdkt, zdk1t )
	342	CALL wrk_dealloc( jpi, jpk, zftw, zdit, zdjt, zdj1t )
[2715]	343	!
[3294]	344	IF( nn_timing == 1 ) CALL timing_stop('ldfght')
	345	!
[3]	346	END SUBROUTINE ldfght
	347
	348	#else
	349	!!----------------------------------------------------------------------
	350	!! Dummy module : NO rotation of the lateral mixing tensor
	351	!!----------------------------------------------------------------------
	352	CONTAINS
[3294]	353	SUBROUTINE tra_ldf_bilapg( kt, kit000, cdtype, ptb, pta, kjpt ) ! Empty routine
	354	INTEGER :: kt, kit000
[2528]	355	CHARACTER(len=3) :: cdtype
	356	REAL, DIMENSION(:,:,:,:) :: ptb, pta
[3294]	357	WRITE(,) 'tra_ldf_iso: You should not have seen this print! error?', &
	358	& kt, kit000, cdtype, ptb(1,1,1,1), pta(1,1,1,1), kjpt
[3]	359	END SUBROUTINE tra_ldf_bilapg
	360	#endif
	361
	362	!!==============================================================================
	363	END MODULE traldf_bilapg

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

Download in other formats: