Context Navigation

traldf_bilapg.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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,*) '~~~~~~~~~~~~~~'
[11101]	90	IF(lwp .AND. lflush) CALL flush(numout)
[3]	91	ENDIF
	92
[2528]	93	! 1. Laplacian of ptb * aht
[3]	94	! -----------------------------
[3294]	95	CALL ldfght( kt, cdtype, ptb, zwk1, kjpt, 1 ) ! rotated harmonic operator applied to ptb and multiply by aht
[2528]	96	! ! output in wk1
	97	!
	98	DO jn = 1, kjpt
[3294]	99	CALL lbc_lnk( zwk1(:,:,:,jn) , 'T', 1. ) ! Lateral boundary conditions on wk1 (unchanged sign)
[2528]	100	END DO
[3]	101
[2528]	102	! 2. Bilaplacian of ptb
[3]	103	! -------------------------
[3294]	104	CALL ldfght( kt, cdtype, zwk1, zwk2, kjpt, 2 ) ! rotated harmonic operator applied to wk1 ; output in wk2
[3]	105
	106
	107	! 3. Update the tracer trends (j-slab : 2, jpj-1)
	108	! ---------------------------
[2528]	109	DO jn = 1, kjpt
	110	DO jj = 2, jpjm1
	111	DO jk = 1, jpkm1
	112	DO ji = 2, jpim1
	113	! add it to the general tracer trends
[3294]	114	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zwk2(ji,jj,jk,jn)
[2528]	115	END DO
[3]	116	END DO
	117	END DO
[2528]	118	END DO
	119	!
[3294]	120	CALL wrk_dealloc( jpi, jpj, jpk, kjpt, zwk1, zwk2 )
[2715]	121	!
[3294]	122	IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_bilapg')
	123	!
[3]	124	END SUBROUTINE tra_ldf_bilapg
	125
	126
[2528]	127	SUBROUTINE ldfght ( kt, cdtype, pt, plt, kjpt, kaht )
[3]	128	!!----------------------------------------------------------------------
	129	!! * ROUTINE ldfght *
	130	!!
[2528]	131	!! ** Purpose : Apply a geopotential harmonic operator to (pt) and
[3]	132	!! multiply it by the eddy diffusivity coefficient (if kaht=1).
	133	!! Routine only used in s-coordinates (l_sco=T) with bilaplacian
	134	!! operator (ln_traldf_bilap=T) acting along geopotential surfaces
	135	!! (ln_traldf_hor).
	136	!!
	137	!! ** Method : The harmonic operator rotated along geopotential
[2528]	138	!! surfaces is applied to (pt) using the slopes of geopotential
[3]	139	!! surfaces computed in inildf routine. The result is provided in
	140	!! (plt,pls) arrays. It is computed in 2 steps:
	141	!!
	142	!! First step: horizontal part of the operator. It is computed on
	143	!! ========== pt as follows (idem on ps)
	144	!! horizontal fluxes :
	145	!! zftu = e2ue3u/e1u di[ pt ] - e2uuslp dk[ mi(mk(pt)) ]
	146	!! zftv = e1ve3v/e2v dj[ pt ] - e1vvslp dk[ mj(mk(pt)) ]
	147	!! take the horizontal divergence of the fluxes (no divided by
	148	!! the volume element :
	149	!! plt = di-1[ zftu ] + dj-1[ zftv ]
	150	!!
	151	!! Second step: vertical part of the operator. It is computed on
	152	!! =========== pt as follows (idem on ps)
	153	!! vertical fluxes :
	154	!! zftw = e1te2t/e3w (wslpi^2+wslpj^2) dk-1[ pt ]
	155	!! - e2t * wslpi di[ mi(mk(pt)) ]
	156	!! - e1t * wslpj dj[ mj(mk(pt)) ]
	157	!! take the vertical divergence of the fluxes add it to the hori-
	158	!! zontal component, divide the result by the volume element and
	159	!! if kaht=1, multiply by the eddy diffusivity coefficient:
	160	!! plt = aht / (e1te2te3t) { plt + dk[ zftw ] }
	161	!! else:
	162	!! plt = 1 / (e1te2te3t) { plt + dk[ zftw ] }
	163	!!
[2528]	164	!!----------------------------------------------------------------------
[2715]	165	USE oce , ONLY: zftv => ua ! ua used as workspace
	166	!
[2528]	167	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	168	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	169	INTEGER , INTENT(in ) :: kjpt !: dimension of
	170	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: pt ! tracer fields ( before for 1st call
	171	! ! and laplacian of these fields for 2nd call.
	172	REAL(wp) , INTENT(out), DIMENSION(jpi,jpj,jpk,kjpt) :: plt !: partial harmonic operator applied to pt components except
	173	! !: second order vertical derivative term
	174	INTEGER , INTENT(in ) :: kaht !: =1 multiply the laplacian by the eddy diffusivity coeff.
	175	! !: =2 no multiplication
	176	!!
	177	INTEGER :: ji, jj, jk,jn ! dummy loop indices
	178	! ! temporary scalars
	179	REAL(wp) :: zabe1, zabe2, zmku, zmkv
	180	REAL(wp) :: zbtr, ztah, ztav
	181	REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4
[3294]	182	REAL(wp), POINTER, DIMENSION(:,:) :: zftu, zdkt, zdk1t
	183	REAL(wp), POINTER, DIMENSION(:,:) :: zftw, zdit, zdjt, zdj1t
[3]	184	!!----------------------------------------------------------------------
[2528]	185	!
[3294]	186	IF( nn_timing == 1 ) CALL timing_start('ldfght')
	187	!
	188	CALL wrk_alloc( jpi, jpj, zftu, zdkt, zdk1t )
	189	CALL wrk_alloc( jpi, jpk, zftw, zdit, zdjt, zdj1t )
	190	!
[2528]	191	DO jn = 1, kjpt
	192	! ! ********** ! ! ===============
	193	DO jk = 1, jpkm1 ! First step ! ! Horizontal slab
	194	! ! ********** ! ! ===============
	195
	196	! I.1 Vertical gradient of pt and ps at level jk and jk+1
	197	! -------------------------------------------------------
	198	! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
	199
	200	zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
	201	IF( jk == 1 ) THEN
	202	zdkt(:,:) = zdk1t(:,:)
	203	ELSE
	204	zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
	205	ENDIF
[3]	206
	207
[2528]	208	! I.2 Horizontal fluxes
	209	! ---------------------
	210
	211	DO jj = 1, jpjm1
	212	DO ji = 1, jpim1
[4292]	213	zabe1 = re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
	214	zabe2 = re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
[2528]	215
	216	zmku = 1./MAX( tmask(ji+1,jj,jk )+tmask(ji,jj,jk+1) &
	217	& +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk ),1. )
	218	zmkv = 1./MAX( tmask(ji,jj+1,jk )+tmask(ji,jj,jk+1) &
	219	& +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk ),1. )
	220
	221	zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
	222	zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
	223
	224	zftu(ji,jj)= umask(ji,jj,jk) * &
	225	& ( zabe1 *( pt (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) ) &
	226	& + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj) &
	227	& +zdk1t(ji+1,jj) + zdkt (ji,jj) ) )
	228
	229	zftv(ji,jj,jk)= vmask(ji,jj,jk) * &
	230	& ( zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) &
	231	& + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj) &
	232	& +zdk1t(ji,jj+1) + zdkt (ji,jj) ) )
	233	END DO
[3]	234	END DO
	235
	236
[2528]	237	! I.3 Second derivative (divergence) (not divided by the volume)
	238	! ---------------------
	239
	240	DO jj = 2 , jpjm1
	241	DO ji = 2 , jpim1
	242	ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
	243	plt(ji,jj,jk,jn) = ztah
	244	END DO
[3]	245	END DO
[2528]	246	! ! ===============
	247	END DO ! End of slab
	248	! ! ===============
	249	! "Poleward" diffusive heat or salt transport
[7179]	250	! note sign is reversed to give down-gradient diffusive transports (#1043)
	251	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) ) CALL dia_ptr_ohst_components( jn, 'ldf', -zftv(:,:,:) )
[3]	252
[2528]	253	! ! ************ ! ! ===============
	254	DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab
	255	! ! ************ ! ! ===============
	256
	257	! II.1 horizontal tracer gradient
	258	! -------------------------------
	259
	260	DO jk = 1, jpk
	261	DO ji = 1, jpim1
	262	zdit (ji,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj ,jk,jn) ) * umask(ji,jj ,jk)
	263	zdjt (ji,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj ,jk,jn) ) * vmask(ji,jj ,jk)
	264	zdj1t(ji,jk) = ( pt(ji ,jj ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
	265	END DO
[3]	266	END DO
	267
	268
[2528]	269	! II.2 Vertical fluxes
	270	! --------------------
	271
	272	! Surface and bottom vertical fluxes set to zero
	273	zftw(:, 1 ) = 0.e0
	274	zftw(:,jpk) = 0.e0
	275
	276	! interior (2=<jk=<jpk-1)
	277	DO jk = 2, jpkm1
	278	DO ji = 2, jpim1
[4292]	279	zcof0 = e12t(ji,jj) / fse3w_n(ji,jj,jk) &
[2528]	280	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	281	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
	282
	283	zmku = 1./MAX( umask(ji ,jj,jk-1)+umask(ji-1,jj,jk) &
	284	& +umask(ji-1,jj,jk-1)+umask(ji ,jj,jk), 1. )
	285
	286	zmkv = 1./MAX( vmask(ji,jj ,jk-1)+vmask(ji,jj-1,jk) &
	287	& +vmask(ji,jj-1,jk-1)+vmask(ji,jj ,jk), 1. )
	288
	289	zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
	290	zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
	291
	292	zftw(ji,jk) = tmask(ji,jj,jk) * &
	293	& ( zcof0 * ( pt (ji,jj,jk-1,jn) - pt (ji ,jj,jk,jn) ) &
	294	& + zcof3 * ( zdit (ji ,jk-1 ) + zdit (ji-1,jk ) &
	295	& +zdit (ji-1 ,jk-1 ) + zdit (ji ,jk ) ) &
	296	& + zcof4 * ( zdjt (ji ,jk-1 ) + zdj1t(ji ,jk) &
	297	& +zdj1t(ji ,jk-1 ) + zdjt (ji ,jk ) ) )
	298	END DO
[3]	299	END DO
	300
	301
[2528]	302	! II.3 Divergence of vertical fluxes added to the horizontal divergence
	303	! ---------------------------------------------------------------------
	304
	305	IF( kaht == 1 ) THEN
	306	! multiply the laplacian by the eddy diffusivity coefficient
	307	DO jk = 1, jpkm1
	308	DO ji = 2, jpim1
	309	! eddy coef. divided by the volume element
[4292]	310	zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
[2528]	311	! vertical divergence
	312	ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
	313	! harmonic operator applied to (pt,ps) and multiply by aht
	314	plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
	315	END DO
[3]	316	END DO
[2528]	317	ELSEIF( kaht == 2 ) THEN
	318	! second call, no multiplication
	319	DO jk = 1, jpkm1
	320	DO ji = 2, jpim1
	321	! inverse of the volume element
[4292]	322	zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
[2528]	323	! vertical divergence
	324	ztav = zftw(ji,jk) - zftw(ji,jk+1)
	325	! harmonic operator applied to (pt,ps)
	326	plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
	327	END DO
[3]	328	END DO
[2528]	329	ELSE
[8280]	330	WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
	331	WRITE(numout,*) ' We stop'
	332	CALL ctl_stop( 'STOP', 'ldfght : unexpected kaht value')
[2528]	333	ENDIF
	334	! ! ===============
	335	END DO ! End of slab
	336	! ! ===============
	337	END DO
	338	!
[3294]	339	CALL wrk_dealloc( jpi, jpj, zftu, zdkt, zdk1t )
	340	CALL wrk_dealloc( jpi, jpk, zftw, zdit, zdjt, zdj1t )
[2715]	341	!
[3294]	342	IF( nn_timing == 1 ) CALL timing_stop('ldfght')
	343	!
[3]	344	END SUBROUTINE ldfght
	345
	346	#else
	347	!!----------------------------------------------------------------------
	348	!! Dummy module : NO rotation of the lateral mixing tensor
	349	!!----------------------------------------------------------------------
	350	CONTAINS
[3294]	351	SUBROUTINE tra_ldf_bilapg( kt, kit000, cdtype, ptb, pta, kjpt ) ! Empty routine
	352	INTEGER :: kt, kit000
[2528]	353	CHARACTER(len=3) :: cdtype
	354	REAL, DIMENSION(:,:,:,:) :: ptb, pta
[3294]	355	WRITE(,) 'tra_ldf_iso: You should not have seen this print! error?', &
	356	& kt, kit000, cdtype, ptb(1,1,1,1), pta(1,1,1,1), kjpt
[3]	357	END SUBROUTINE tra_ldf_bilapg
	358	#endif
	359
	360	!!==============================================================================
	361	END MODULE traldf_bilapg

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 11101

Download in other formats: