Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

traldf_bilapg.F90 @ 9987

Last change on this file since 9987 was 9987, checked in by emmafiedler, 6 years ago
Merge with GO6 FOAMv14 package branch r9288
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
[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
[9987]	249	! note sign is reversed to give down-gradient diffusive transports (#1043)
	250	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) ) CALL dia_ptr_ohst_components( jn, 'ldf', -zftv(:,:,:) )
[3]	251
[2528]	252	! ! ************ ! ! ===============
	253	DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab
	254	! ! ************ ! ! ===============
	255
	256	! II.1 horizontal tracer gradient
	257	! -------------------------------
	258
	259	DO jk = 1, jpk
	260	DO ji = 1, jpim1
	261	zdit (ji,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj ,jk,jn) ) * umask(ji,jj ,jk)
	262	zdjt (ji,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj ,jk,jn) ) * vmask(ji,jj ,jk)
	263	zdj1t(ji,jk) = ( pt(ji ,jj ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
	264	END DO
[3]	265	END DO
	266
	267
[2528]	268	! II.2 Vertical fluxes
	269	! --------------------
	270
	271	! Surface and bottom vertical fluxes set to zero
	272	zftw(:, 1 ) = 0.e0
	273	zftw(:,jpk) = 0.e0
	274
	275	! interior (2=<jk=<jpk-1)
	276	DO jk = 2, jpkm1
	277	DO ji = 2, jpim1
[4292]	278	zcof0 = e12t(ji,jj) / fse3w_n(ji,jj,jk) &
[2528]	279	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	280	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
	281
	282	zmku = 1./MAX( umask(ji ,jj,jk-1)+umask(ji-1,jj,jk) &
	283	& +umask(ji-1,jj,jk-1)+umask(ji ,jj,jk), 1. )
	284
	285	zmkv = 1./MAX( vmask(ji,jj ,jk-1)+vmask(ji,jj-1,jk) &
	286	& +vmask(ji,jj-1,jk-1)+vmask(ji,jj ,jk), 1. )
	287
	288	zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
	289	zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
	290
	291	zftw(ji,jk) = tmask(ji,jj,jk) * &
	292	& ( zcof0 * ( pt (ji,jj,jk-1,jn) - pt (ji ,jj,jk,jn) ) &
	293	& + zcof3 * ( zdit (ji ,jk-1 ) + zdit (ji-1,jk ) &
	294	& +zdit (ji-1 ,jk-1 ) + zdit (ji ,jk ) ) &
	295	& + zcof4 * ( zdjt (ji ,jk-1 ) + zdj1t(ji ,jk) &
	296	& +zdj1t(ji ,jk-1 ) + zdjt (ji ,jk ) ) )
	297	END DO
[3]	298	END DO
	299
	300
[2528]	301	! II.3 Divergence of vertical fluxes added to the horizontal divergence
	302	! ---------------------------------------------------------------------
	303
	304	IF( kaht == 1 ) THEN
	305	! multiply the laplacian by the eddy diffusivity coefficient
	306	DO jk = 1, jpkm1
	307	DO ji = 2, jpim1
	308	! eddy coef. divided by the volume element
[4292]	309	zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
[2528]	310	! vertical divergence
	311	ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
	312	! harmonic operator applied to (pt,ps) and multiply by aht
	313	plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
	314	END DO
[3]	315	END DO
[2528]	316	ELSEIF( kaht == 2 ) THEN
	317	! second call, no multiplication
	318	DO jk = 1, jpkm1
	319	DO ji = 2, jpim1
	320	! inverse of the volume element
[4292]	321	zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
[2528]	322	! vertical divergence
	323	ztav = zftw(ji,jk) - zftw(ji,jk+1)
	324	! harmonic operator applied to (pt,ps)
	325	plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
	326	END DO
[3]	327	END DO
[2528]	328	ELSE
[9987]	329	WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
	330	WRITE(numout,*) ' We stop'
	331	CALL ctl_stop( 'STOP', 'ldfght : unexpected kaht value')
[2528]	332	ENDIF
	333	! ! ===============
	334	END DO ! End of slab
	335	! ! ===============
	336	END DO
	337	!
[3294]	338	CALL wrk_dealloc( jpi, jpj, zftu, zdkt, zdk1t )
	339	CALL wrk_dealloc( jpi, jpk, zftw, zdit, zdjt, zdj1t )
[2715]	340	!
[3294]	341	IF( nn_timing == 1 ) CALL timing_stop('ldfght')
	342	!
[3]	343	END SUBROUTINE ldfght
	344
	345	#else
	346	!!----------------------------------------------------------------------
	347	!! Dummy module : NO rotation of the lateral mixing tensor
	348	!!----------------------------------------------------------------------
	349	CONTAINS
[3294]	350	SUBROUTINE tra_ldf_bilapg( kt, kit000, cdtype, ptb, pta, kjpt ) ! Empty routine
	351	INTEGER :: kt, kit000
[2528]	352	CHARACTER(len=3) :: cdtype
	353	REAL, DIMENSION(:,:,:,:) :: ptb, pta
[3294]	354	WRITE(,) 'tra_ldf_iso: You should not have seen this print! error?', &
	355	& kt, kit000, cdtype, ptb(1,1,1,1), pta(1,1,1,1), kjpt
[3]	356	END SUBROUTINE tra_ldf_bilapg
	357	#endif
	358
	359	!!==============================================================================
	360	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_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 9987

Download in other formats: