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.

trdicp.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago
Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
Property svn:keywords set to `Id`
File size: 34.3 KB

Rev	Line
[215]	1	MODULE trdicp
	2	!!======================================================================
	3	!! * MODULE trdicp *
	4	!! Ocean diagnostics: ocean tracers and dynamic trends
	5	!!=====================================================================
[2528]	6	!! History : 1.0 ! 2004-08 (C. Talandier) New trends organization
[503]	7	!!----------------------------------------------------------------------
[215]	8	#if defined key_trdtra \|\| defined key_trddyn \|\| defined key_esopa
	9	!!----------------------------------------------------------------------
	10	!! 'key_trdtra' or active tracers trends diagnostics
	11	!! 'key_trddyn' momentum trends diagnostics
	12	!!----------------------------------------------------------------------
[503]	13	!! trd_icp : compute the basin averaged properties for tra/dyn
[215]	14	!! trd_dwr : print dynmaic trends in ocean.output file
	15	!! trd_twr : print tracers trends in ocean.output file
	16	!! trd_icp_init : initialization step
	17	!!----------------------------------------------------------------------
	18	USE oce ! ocean dynamics and tracers variables
	19	USE dom_oce ! ocean space and time domain variables
	20	USE trdmod_oce ! ocean variables trends
	21	USE ldftra_oce ! ocean active tracers: lateral physics
	22	USE ldfdyn_oce ! ocean dynamics: lateral physics
	23	USE zdf_oce ! ocean vertical physics
	24	USE in_out_manager ! I/O manager
	25	USE lib_mpp ! distibuted memory computing library
	26	USE eosbn2 ! equation of state
	27	USE phycst ! physical constants
	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[503]	32	INTERFACE trd_icp
[215]	33	MODULE PROCEDURE trd_2d, trd_3d
	34	END INTERFACE
	35
[503]	36	PUBLIC trd_icp ! called by trdmod.F90
	37	PUBLIC trd_dwr ! called by step.F90
	38	PUBLIC trd_twr ! called by step.F90
	39	PUBLIC trd_icp_init ! called by opa.F90
[215]	40
[3211]	41	!! * Control permutation of array indices
	42	# include "oce_ftrans.h90"
	43	# include "dom_oce_ftrans.h90"
	44	# include "trdmld_oce_ftrans.h90"
	45	# include "ldftra_oce_ftrans.h90"
	46	# include "ldfdyn_oce_ftrans.h90"
	47	# include "zdf_oce_ftrans.h90"
	48
[215]	49	!! * Substitutions
	50	# include "domzgr_substitute.h90"
	51	# include "vectopt_loop_substitute.h90"
	52	!!----------------------------------------------------------------------
[2528]	53	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	54	!! $Id$
[2715]	55	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[215]	56	!!----------------------------------------------------------------------
	57	CONTAINS
	58
[2528]	59	SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )
[215]	60	!!---------------------------------------------------------------------
	61	!! * ROUTINE trd_2d *
	62	!!
	63	!! ** Purpose : verify the basin averaged properties of tracers and/or
[1601]	64	!! momentum equations at every time step frequency nn_trd.
[503]	65	!!----------------------------------------------------------------------
[2715]	66	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptrd2dx ! Temperature or U trend
	67	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptrd2dy ! Salinity or V trend
	68	INTEGER , INTENT(in ) :: ktrd ! tracer trend index
	69	CHARACTER(len=3) , INTENT(in ) :: ctype ! momentum ('DYN') or tracers ('TRA') trends
[215]	70	!!
[2715]	71	INTEGER :: ji, jj ! loop indices
[215]	72	!!----------------------------------------------------------------------
	73
[2715]	74	SELECT CASE( ctype ) !== Mask trends ==!
[503]	75	!
[2715]	76	CASE( 'DYN' ) ! Momentum
[215]	77	DO jj = 1, jpjm1
	78	DO ji = 1, jpim1
[2715]	79	ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * tmask_i(ji+1,jj ) * tmask_i(ji,jj) * umask(ji,jj,1)
	80	ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * tmask_i(ji ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
[215]	81	END DO
	82	END DO
[2715]	83	ptrd2dx(jpi, : ) = 0._wp ; ptrd2dy(jpi, : ) = 0._wp
	84	ptrd2dx( : ,jpj) = 0._wp ; ptrd2dy( : ,jpj) = 0._wp
[503]	85	!
[2715]	86	CASE( 'TRA' ) ! Tracers
[215]	87	ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
	88	ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
[503]	89	!
[215]	90	END SELECT
	91
[2715]	92	SELECT CASE( ctype ) !== Basin averaged tracer/momentum trends ==!
[503]	93	!
[2715]	94	CASE( 'DYN' ) ! Momentum
	95	umo(ktrd) = 0._wp
	96	vmo(ktrd) = 0._wp
[503]	97	!
	98	SELECT CASE( ktrd )
	99	CASE( jpdyn_trd_swf ) ! surface forcing
[2715]	100	umo(ktrd) = SUM( ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1) )
	101	vmo(ktrd) = SUM( ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1) )
[215]	102	END SELECT
[503]	103	!
	104	CASE( 'TRA' ) ! Tracers
[2715]	105	tmo(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
	106	smo(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
[215]	107	END SELECT
	108
[2715]	109	SELECT CASE( ctype ) !== Basin averaged tracer/momentum square trends ==! (now field)
[503]	110	!
	111	CASE( 'DYN' ) ! Momentum
[2715]	112	hke(ktrd) = SUM( un(:,:,1) * ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1) &
	113	& + vn(:,:,1) * ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1) )
[503]	114	!
	115	CASE( 'TRA' ) ! Tracers
[2715]	116	t2(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tn(:,:,1) )
	117	s2(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) * sn(:,:,1) )
[503]	118	!
[215]	119	END SELECT
[503]	120	!
[215]	121	END SUBROUTINE trd_2d
	122
	123
[2528]	124	SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd, ctype )
[215]	125	!!---------------------------------------------------------------------
	126	!! * ROUTINE trd_3d *
	127	!!
	128	!! ** Purpose : verify the basin averaged properties of tracers and/or
[1601]	129	!! momentum equations at every time step frequency nn_trd.
[503]	130	!!----------------------------------------------------------------------
[3211]	131
	132	!! DCSE_NEMO: This style defeats ftrans
	133	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptrd3dx ! Temperature or U trend
	134	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptrd3dy ! Salinity or V trend
	135
	136	!FTRANS ptrd3dx ptrd3dy :I :I :z
[4409]	137	REAL(wp), INTENT(inout) :: ptrd3dx(jpi,jpj,jpkorig) ! Temperature or U trend
	138	REAL(wp), INTENT(inout) :: ptrd3dy(jpi,jpj,jpkorig) ! Salinity or V trend
[3211]	139
[2715]	140	INTEGER, INTENT(in ) :: ktrd ! momentum or tracer trend index
	141	CHARACTER(len=3), INTENT(in ) :: ctype ! momentum ('DYN') or tracers ('TRA') trends
[215]	142	!!
[2715]	143	INTEGER :: ji, jj, jk ! dummy loop indices
[215]	144	!!----------------------------------------------------------------------
	145
[2715]	146	SELECT CASE( ctype ) !== Mask the trends ==!
[503]	147	!
	148	CASE( 'DYN' ) ! Momentum
[3211]	149	#if defined key_z_first
	150	DO jj = 1, jpjm1
	151	DO ji = 1, jpim1
	152	DO jk = 1, jpkm1
	153	#else
[2715]	154	DO jk = 1, jpkm1
[215]	155	DO jj = 1, jpjm1
	156	DO ji = 1, jpim1
[3211]	157	#endif
[2715]	158	ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * tmask_i(ji+1,jj ) * tmask_i(ji,jj) * umask(ji,jj,jk)
	159	ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * tmask_i(ji ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
[503]	160	END DO
	161	END DO
	162	END DO
[2715]	163	ptrd3dx(jpi, : ,:) = 0._wp ; ptrd3dy(jpi, : ,:) = 0._wp
	164	ptrd3dx( : ,jpj,:) = 0._wp ; ptrd3dy( : ,jpj,:) = 0._wp
[503]	165	!
	166	CASE( 'TRA' ) ! Tracers
[3211]	167	#if defined key_z_first
	168	DO jj = 1, jpj
	169	DO ji = 1, jpi
	170	DO jk = 1, jpkm1
	171	ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
	172	ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
	173	END DO
	174	END DO
	175	END DO
	176	#else
[2715]	177	DO jk = 1, jpkm1
[215]	178	ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
	179	ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
[503]	180	END DO
[3211]	181	#endif
[503]	182	!
[215]	183	END SELECT
	184
[2715]	185	SELECT CASE( ctype ) !== Basin averaged tracer/momentum trends ==!
[503]	186	!
	187	CASE( 'DYN' ) ! Momentum
[2715]	188	umo(ktrd) = 0._wp
	189	vmo(ktrd) = 0._wp
[3211]	190	#if defined key_z_first
	191	!! DCSE_NEMO: this changes the order of summation
	192	DO jj = 1, jpj
	193	DO ji = 1, jpi
	194	DO jk = 1, jpkm1
	195	umo(ktrd) = umo(ktrd) + ptrd3dx(ji,jj,jk) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
	196	vmo(ktrd) = vmo(ktrd) + ptrd3dy(ji,jj,jk) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	197	END DO
	198	END DO
	199	END DO
	200	#else
[2715]	201	DO jk = 1, jpkm1
	202	umo(ktrd) = umo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk) )
	203	vmo(ktrd) = vmo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk) )
[215]	204	END DO
[3211]	205	#endif
[503]	206	!
	207	CASE( 'TRA' ) ! Tracers
[2715]	208	tmo(ktrd) = 0._wp
	209	smo(ktrd) = 0._wp
[3211]	210	#if defined key_z_first
	211	!! DCSE_NEMO: this changes the order of summation
	212	DO jj = 1, jpj
	213	DO ji = 1, jpi
	214	DO jk = 1, jpkm1
	215	tmo(ktrd) = tmo(ktrd) + ptrd3dx(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk) )
	216	smo(ktrd) = smo(ktrd) + ptrd3dy(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk) )
	217	END DO
	218	END DO
	219	END DO
	220	#else
[215]	221	DO jk = 1, jpkm1
[2715]	222	tmo(ktrd) = tmo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
	223	smo(ktrd) = smo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
[215]	224	END DO
[3211]	225	#endif
[503]	226	!
[215]	227	END SELECT
	228
[2715]	229	SELECT CASE( ctype ) !== Basin averaged tracer/momentum square trends ==! (now field)
[503]	230	!
	231	CASE( 'DYN' ) ! Momentum
[2715]	232	hke(ktrd) = 0._wp
[3211]	233	#if defined key_z_first
	234	!! DCSE_NEMO: this changes the order of summation
	235	DO jj = 1, jpj
	236	DO ji = 1, jpi
	237	DO jk = 1, jpkm1
	238	hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * ptrd3dx(ji,jj,jk) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
	239	& + vn(ji,jj,jk) * ptrd3dy(ji,jj,jk) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	240	END DO
	241	END DO
	242	END DO
	243	#else
[2715]	244	DO jk = 1, jpkm1
	245	hke(ktrd) = hke(ktrd) + SUM( un(:,:,jk) * ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk) &
	246	& + vn(:,:,jk) * ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk) )
[215]	247	END DO
[3211]	248	#endif
[503]	249	!
	250	CASE( 'TRA' ) ! Tracers
[2715]	251	t2(ktrd) = 0._wp
	252	s2(ktrd) = 0._wp
[3211]	253	#if defined key_z_first
	254	!! DCSE_NEMO: this changes the order of summation
	255	DO jj = 1, jpj
	256	DO ji = 1, jpi
	257	DO jk = 1, jpkm1
	258	t2(ktrd) = t2(ktrd) + ptrd3dx(ji,jj,jk) * tn(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk)
	259	s2(ktrd) = s2(ktrd) + ptrd3dy(ji,jj,jk) * sn(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk)
	260	END DO
	261	END DO
	262	END DO
	263	#else
[2715]	264	DO jk = 1, jpkm1
[3211]	265	t2(ktrd) = t2(ktrd) + SUM( ptrd3dx(:,:,jk) * tn(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
	266	s2(ktrd) = s2(ktrd) + SUM( ptrd3dy(:,:,jk) * sn(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
[215]	267	END DO
[3211]	268	#endif
[503]	269	!
[215]	270	END SELECT
[503]	271	!
[215]	272	END SUBROUTINE trd_3d
	273
	274
	275	SUBROUTINE trd_icp_init
	276	!!---------------------------------------------------------------------
	277	!! * ROUTINE trd_icp_init *
	278	!!
[503]	279	!! ** Purpose : Read the namtrd namelist
[215]	280	!!----------------------------------------------------------------------
[2715]	281	INTEGER :: ji, jj, jk ! dummy loop indices
[215]	282	!!----------------------------------------------------------------------
	283
	284	IF(lwp) THEN
	285	WRITE(numout,*)
	286	WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
	287	WRITE(numout,*) '~~~~~~~~~~~~~'
	288	ENDIF
	289
	290	! Total volume at t-points:
[2715]	291	tvolt = 0._wp
[3211]	292	#if defined key_z_first
	293	DO jj = 1, jpj
	294	DO ji = 1, jpi
	295	DO jk = 1, jpkm1
	296	tvolt = tvolt + e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
	297	END DO
	298	END DO
	299	END DO
	300	#else
[215]	301	DO jk = 1, jpkm1
[3211]	302	!! DCSE_NEMO: This looks plain wrong
	303	! tvolt = SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
	304	tvolt = tvolt + SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
[215]	305	END DO
[3211]	306	#endif
[215]	307	IF( lk_mpp ) CALL mpp_sum( tvolt ) ! sum over the global domain
	308
[503]	309	IF(lwp) WRITE(numout,*) ' total ocean volume at T-point tvolt = ',tvolt
[215]	310
	311	#if defined key_trddyn
	312	! Initialization of potential to kinetic energy conversion
[2715]	313	rpktrd = 0._wp
[215]	314
	315	! Total volume at u-, v- points:
[2715]	316	tvolu = 0._wp
	317	tvolv = 0._wp
[215]	318
[3211]	319	#if defined key_z_first
	320	DO jj = 2, jpjm1
	321	DO ji = 2, jpim1
	322	DO jk = 1, jpk
	323	#else
[215]	324	DO jk = 1, jpk
	325	DO jj = 2, jpjm1
	326	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	327	#endif
[2715]	328	tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) * tmask_i(ji+1,jj ) * tmask_i(ji,jj) * umask(ji,jj,jk)
	329	tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
[215]	330	END DO
	331	END DO
	332	END DO
	333	IF( lk_mpp ) CALL mpp_sum( tvolu ) ! sums over the global domain
	334	IF( lk_mpp ) CALL mpp_sum( tvolv )
	335
	336	IF(lwp) THEN
[503]	337	WRITE(numout,*) ' total ocean volume at U-point tvolu = ',tvolu
	338	WRITE(numout,*) ' total ocean volume at V-point tvolv = ',tvolv
[215]	339	ENDIF
	340	#endif
[503]	341	!
[215]	342	END SUBROUTINE trd_icp_init
	343
	344
	345	SUBROUTINE trd_dwr( kt )
	346	!!---------------------------------------------------------------------
	347	!! * ROUTINE trd_dwr *
	348	!!
	349	!! ** Purpose : write dynamic trends in ocean.output
[503]	350	!!----------------------------------------------------------------------
[2715]	351	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	352	USE wrk_nemo, ONLY: zkepe => wrk_3d_1 , zkx => wrk_3d_2 ! 3D workspace
	353	USE wrk_nemo, ONLY: zky => wrk_3d_3 , zkz => wrk_3d_4 ! - -
[3211]	354
	355	!! DCSE_NEMO: need additional directives for renamed module variables
	356	!FTRANS zkepe zkx zky zkz :I :I :z
	357
[2715]	358	!
	359	INTEGER, INTENT(in) :: kt ! ocean time-step index
	360	!
	361	INTEGER :: ji, jj, jk ! dummy loop indices
	362	REAL(wp) :: zcof ! local scalar
[215]	363	!!----------------------------------------------------------------------
	364
[2715]	365	IF( wrk_in_use(3, 1,2,3,4) ) THEN
	366	CALL ctl_stop('trd_dwr: requested workspace arrays unavailable') ; RETURN
	367	ENDIF
	368
[215]	369	! I. Momentum trends
	370	! -------------------
	371
[1601]	372	IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
[215]	373
	374	! I.1 Conversion potential energy - kinetic energy
	375	! --------------------------------------------------
	376	! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
[2715]	377	zkx (:,:,:) = 0._wp
	378	zky (:,:,:) = 0._wp
	379	zkz (:,:,:) = 0._wp
	380	zkepe(:,:,:) = 0._wp
[215]	381
[2528]	382	CALL eos( tsn, rhd, rhop ) ! now potential and in situ densities
[215]	383
[2715]	384	zcof = 0.5_wp / rau0 ! Density flux at w-point
[3211]	385	#if defined key_z_first
	386	DO jj = 1, jpj
	387	DO ji = 1, jpi
	388	zkz(ji,jj,1) = 0._wp
	389	DO jk = 2, jpk
	390	zkz(ji,jj,jk) = e1e2t(ji,jj) * wn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) ) * tmask_i(ji,jj)
	391	END DO
	392	END DO
	393	END DO
	394	#else
[2715]	395	zkz(:,:,1) = 0._wp
[215]	396	DO jk = 2, jpk
[2715]	397	zkz(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
[215]	398	END DO
[3211]	399	#endif
[215]	400
[2715]	401	zcof = 0.5_wp / rau0 ! Density flux at u and v-points
[3211]	402	#if defined key_z_first
	403	DO jj = 1, jpjm1
	404	DO ji = 1, jpim1
	405	DO jk = 1, jpkm1
	406	#else
[2715]	407	DO jk = 1, jpkm1
[215]	408	DO jj = 1, jpjm1
	409	DO ji = 1, jpim1
[3211]	410	#endif
[2715]	411	zkx(ji,jj,jk) = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
	412	zky(ji,jj,jk) = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
[215]	413	END DO
	414	END DO
	415	END DO
	416
[3211]	417	#if defined key_z_first
	418	DO jj = 2, jpjm1 ! Density flux divergence at t-point
	419	DO ji = 2, jpim1
	420	DO jk = 1, jpkm1
	421	#else
[2715]	422	DO jk = 1, jpkm1 ! Density flux divergence at t-point
[215]	423	DO jj = 2, jpjm1
	424	DO ji = 2, jpim1
[3211]	425	#endif
[2715]	426	zkepe(ji,jj,jk) = - ( zkz(ji,jj,jk) - zkz(ji ,jj ,jk+1) &
	427	& + zkx(ji,jj,jk) - zkx(ji-1,jj ,jk ) &
	428	& + zky(ji,jj,jk) - zky(ji ,jj-1,jk ) ) &
	429	& / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
[215]	430	END DO
	431	END DO
	432	END DO
	433
	434	! I.2 Basin averaged kinetic energy trend
	435	! ----------------------------------------
[2715]	436	peke = 0._wp
[3211]	437	#if defined key_z_first
	438	DO jj = 1, jpj
	439	DO ji = 1, jpi
	440	DO jk = 1, jpkm1
	441	peke = peke + zkepe(ji,jj,jk) * fsdept(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk)
	442	END DO
	443	END DO
	444	END DO
	445	#else
[2715]	446	DO jk = 1, jpkm1
	447	peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
[215]	448	END DO
[3211]	449	#endif
[2715]	450	peke = grav * peke
[215]	451
	452	! I.3 Sums over the global domain
	453	! ---------------------------------
	454	IF( lk_mpp ) THEN
[503]	455	CALL mpp_sum( peke )
	456	CALL mpp_sum( umo , jptot_dyn )
	457	CALL mpp_sum( vmo , jptot_dyn )
	458	CALL mpp_sum( hke , jptot_dyn )
	459	ENDIF
[215]	460
	461	! I.2 Print dynamic trends in the ocean.output file
	462	! --------------------------------------------------
	463
	464	IF(lwp) THEN
	465	WRITE (numout,*)
	466	WRITE (numout,*)
	467	WRITE (numout,9500) kt
[503]	468	WRITE (numout,9501) umo(jpicpd_hpg) / tvolu, vmo(jpicpd_hpg) / tvolv
	469	WRITE (numout,9502) umo(jpicpd_keg) / tvolu, vmo(jpicpd_keg) / tvolv
	470	WRITE (numout,9503) umo(jpicpd_rvo) / tvolu, vmo(jpicpd_rvo) / tvolv
	471	WRITE (numout,9504) umo(jpicpd_pvo) / tvolu, vmo(jpicpd_pvo) / tvolv
	472	WRITE (numout,9505) umo(jpicpd_ldf) / tvolu, vmo(jpicpd_ldf) / tvolv
[1129]	473	WRITE (numout,9506) umo(jpicpd_had) / tvolu, vmo(jpicpd_had) / tvolv
	474	WRITE (numout,9507) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
	475	WRITE (numout,9508) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
	476	WRITE (numout,9509) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
	477	WRITE (numout,9510) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
	478	WRITE (numout,9511) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
	479	WRITE (numout,9512) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
	480	WRITE (numout,9513)
	481	WRITE (numout,9514) &
[503]	482	& ( umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf) &
[1129]	483	& + umo(jpicpd_had) + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat) &
	484	& + umo(jpicpd_swf) + umo(jpicpd_bfr) ) / tvolu, &
[503]	485	& ( vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf) &
[1129]	486	& + vmo(jpicpd_had) + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat) &
	487	& + vmo(jpicpd_swf) + vmo(jpicpd_bfr) ) / tvolv
[215]	488	ENDIF
	489
	490	9500 FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
	491	9501 FORMAT(' pressure gradient u= ', e20.13, ' v= ', e20.13)
	492	9502 FORMAT(' ke gradient u= ', e20.13, ' v= ', e20.13)
	493	9503 FORMAT(' relative vorticity term u= ', e20.13, ' v= ', e20.13)
	494	9504 FORMAT(' coriolis term u= ', e20.13, ' v= ', e20.13)
	495	9505 FORMAT(' horizontal diffusion u= ', e20.13, ' v= ', e20.13)
[1129]	496	9506 FORMAT(' horizontal advection u= ', e20.13, ' v= ', e20.13)
	497	9507 FORMAT(' vertical advection u= ', e20.13, ' v= ', e20.13)
	498	9508 FORMAT(' vertical diffusion u= ', e20.13, ' v= ', e20.13)
	499	9509 FORMAT(' surface pressure gradient u= ', e20.13, ' v= ', e20.13)
	500	9510 FORMAT(' surface wind forcing u= ', e20.13, ' v= ', e20.13)
	501	9511 FORMAT(' dampimg term u= ', e20.13, ' v= ', e20.13)
	502	9512 FORMAT(' bottom flux u= ', e20.13, ' v= ', e20.13)
	503	9513 FORMAT(' -----------------------------------------------------------------------------')
	504	9514 FORMAT(' total trend u= ', e20.13, ' v= ', e20.13)
[215]	505
	506	IF(lwp) THEN
	507	WRITE (numout,*)
	508	WRITE (numout,*)
	509	WRITE (numout,9520) kt
[503]	510	WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
	511	WRITE (numout,9522) hke(jpicpd_keg) / tvolt
	512	WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
	513	WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
	514	WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
[1129]	515	WRITE (numout,9526) hke(jpicpd_had) / tvolt
	516	WRITE (numout,9527) hke(jpicpd_zad) / tvolt
	517	WRITE (numout,9528) hke(jpicpd_zdf) / tvolt
	518	WRITE (numout,9529) hke(jpicpd_spg) / tvolt
	519	WRITE (numout,9530) hke(jpicpd_swf) / tvolt
	520	WRITE (numout,9531) hke(jpicpd_dat) / tvolt
[1708]	521	WRITE (numout,9532) hke(jpicpd_bfr) / tvolt
	522	WRITE (numout,9533)
	523	WRITE (numout,9534) &
[503]	524	& ( hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf) &
[1129]	525	& + hke(jpicpd_had) + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat) &
[1708]	526	& + hke(jpicpd_swf) + hke(jpicpd_bfr) ) / tvolt
[215]	527	ENDIF
	528
	529	9520 FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
	530	9521 FORMAT(' pressure gradient u2= ', e20.13)
	531	9522 FORMAT(' ke gradient u2= ', e20.13)
	532	9523 FORMAT(' relative vorticity term u2= ', e20.13)
	533	9524 FORMAT(' coriolis term u2= ', e20.13)
	534	9525 FORMAT(' horizontal diffusion u2= ', e20.13)
[1129]	535	9526 FORMAT(' horizontal advection u2= ', e20.13)
	536	9527 FORMAT(' vertical advection u2= ', e20.13)
	537	9528 FORMAT(' vertical diffusion u2= ', e20.13)
	538	9529 FORMAT(' surface pressure gradient u2= ', e20.13)
	539	9530 FORMAT(' surface wind forcing u2= ', e20.13)
	540	9531 FORMAT(' dampimg term u2= ', e20.13)
[1708]	541	9532 FORMAT(' bottom flux u2= ', e20.13)
	542	9533 FORMAT(' --------------------------------------------------')
	543	9534 FORMAT(' total trend u2= ', e20.13)
[215]	544
	545	IF(lwp) THEN
	546	WRITE (numout,*)
	547	WRITE (numout,*)
	548	WRITE (numout,9540) kt
[1129]	549	WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
	550	WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
[503]	551	WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
	552	WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
	553	WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
	554	WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
	555	WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
	556	WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
	557	WRITE (numout,*)
	558	WRITE (numout,*)
[215]	559	ENDIF
	560
	561	9540 FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
	562	9541 FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
[1129]	563	9542 FORMAT(' 0 = ke gradient + horizontal + vertical advection : ', e20.13)
[215]	564	9543 FORMAT(' 0 = coriolis term (true if key_vorenergy or key_combined): ', e20.13)
[503]	565	9544 FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.) : ', e20.13)
	566	9545 FORMAT(' 0 = surface pressure gradient : ', e20.13)
	567	9546 FORMAT(' 0 > horizontal diffusion : ', e20.13)
	568	9547 FORMAT(' 0 > vertical diffusion : ', e20.13)
	569	9548 FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop) : ', e20.13, ' u.dz(rhop) =', e20.13)
	570	!
[215]	571	! Save potential to kinetic energy conversion for next time step
	572	rpktrd = peke
[503]	573	!
[215]	574	ENDIF
[503]	575	!
[2715]	576	IF( wrk_not_released(3, 1,2,3,4) ) CALL ctl_stop('trd_dwr: failed to release workspace arrays')
	577	!
[215]	578	END SUBROUTINE trd_dwr
	579
	580
	581	SUBROUTINE trd_twr( kt )
	582	!!---------------------------------------------------------------------
	583	!! * ROUTINE trd_twr *
	584	!!
	585	!! ** Purpose : write active tracers trends in ocean.output
	586	!!----------------------------------------------------------------------
[503]	587	INTEGER, INTENT(in) :: kt ! ocean time-step index
[215]	588	!!----------------------------------------------------------------------
	589
	590	! I. Tracers trends
	591	! -----------------
	592
[1601]	593	IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
[215]	594
	595	! I.1 Sums over the global domain
	596	! -------------------------------
	597	IF( lk_mpp ) THEN
[503]	598	CALL mpp_sum( tmo, jptot_tra )
	599	CALL mpp_sum( smo, jptot_tra )
	600	CALL mpp_sum( t2 , jptot_tra )
	601	CALL mpp_sum( s2 , jptot_tra )
[215]	602	ENDIF
	603
	604	! I.2 Print tracers trends in the ocean.output file
	605	! --------------------------------------------------
	606
	607	IF(lwp) THEN
	608	WRITE (numout,*)
	609	WRITE (numout,*)
	610	WRITE (numout,9400) kt
[2528]	611	WRITE (numout,9401) tmo(jpicpt_xad) / tvolt, smo(jpicpt_xad) / tvolt
	612	WRITE (numout,9411) tmo(jpicpt_yad) / tvolt, smo(jpicpt_yad) / tvolt
[503]	613	WRITE (numout,9402) tmo(jpicpt_zad) / tvolt, smo(jpicpt_zad) / tvolt
	614	WRITE (numout,9403) tmo(jpicpt_ldf) / tvolt, smo(jpicpt_ldf) / tvolt
	615	WRITE (numout,9404) tmo(jpicpt_zdf) / tvolt, smo(jpicpt_zdf) / tvolt
	616	WRITE (numout,9405) tmo(jpicpt_npc) / tvolt, smo(jpicpt_npc) / tvolt
	617	WRITE (numout,9406) tmo(jpicpt_dmp) / tvolt, smo(jpicpt_dmp) / tvolt
	618	WRITE (numout,9407) tmo(jpicpt_qsr) / tvolt
	619	WRITE (numout,9408) tmo(jpicpt_nsr) / tvolt, smo(jpicpt_nsr) / tvolt
	620	WRITE (numout,9409)
	621	WRITE (numout,9410) ( tmo(jpicpt_xad) + tmo(jpicpt_yad) + tmo(jpicpt_zad) + tmo(jpicpt_ldf) + tmo(jpicpt_zdf) &
	622	& + tmo(jpicpt_npc) + tmo(jpicpt_dmp) + tmo(jpicpt_qsr) + tmo(jpicpt_nsr) ) / tvolt, &
	623	& ( smo(jpicpt_xad) + smo(jpicpt_yad) + smo(jpicpt_zad) + smo(jpicpt_ldf) + smo(jpicpt_zdf) &
	624	& + smo(jpicpt_npc) + smo(jpicpt_dmp) + smo(jpicpt_nsr) ) / tvolt
[215]	625	ENDIF
	626
	627	9400 FORMAT(' tracer trend at it= ',i6,' : temperature', &
	628	' salinity',/' ============================')
[2528]	629	9401 FORMAT(' zonal advection ',e20.13,' ',e20.13)
	630	9411 FORMAT(' meridional advection ',e20.13,' ',e20.13)
[215]	631	9402 FORMAT(' vertical advection ',e20.13,' ',e20.13)
	632	9403 FORMAT(' horizontal diffusion ',e20.13,' ',e20.13)
	633	9404 FORMAT(' vertical diffusion ',e20.13,' ',e20.13)
[503]	634	9405 FORMAT(' static instability mixing ',e20.13,' ',e20.13)
[215]	635	9406 FORMAT(' damping term ',e20.13,' ',e20.13)
[503]	636	9407 FORMAT(' penetrative qsr ',e20.13)
	637	9408 FORMAT(' non solar radiation ',e20.13,' ',e20.13)
[215]	638	9409 FORMAT(' -------------------------------------------------------------------------')
	639	9410 FORMAT(' total trend ',e20.13,' ',e20.13)
	640
	641
	642	IF(lwp) THEN
	643	WRITE (numout,*)
	644	WRITE (numout,*)
	645	WRITE (numout,9420) kt
[2528]	646	WRITE (numout,9421) t2(jpicpt_xad) / tvolt, s2(jpicpt_xad) / tvolt
	647	WRITE (numout,9431) t2(jpicpt_yad) / tvolt, s2(jpicpt_yad) / tvolt
[503]	648	WRITE (numout,9422) t2(jpicpt_zad) / tvolt, s2(jpicpt_zad) / tvolt
	649	WRITE (numout,9423) t2(jpicpt_ldf) / tvolt, s2(jpicpt_ldf) / tvolt
	650	WRITE (numout,9424) t2(jpicpt_zdf) / tvolt, s2(jpicpt_zdf) / tvolt
	651	WRITE (numout,9425) t2(jpicpt_npc) / tvolt, s2(jpicpt_npc) / tvolt
	652	WRITE (numout,9426) t2(jpicpt_dmp) / tvolt, s2(jpicpt_dmp) / tvolt
	653	WRITE (numout,9427) t2(jpicpt_qsr) / tvolt
	654	WRITE (numout,9428) t2(jpicpt_nsr) / tvolt, s2(jpicpt_nsr) / tvolt
[215]	655	WRITE (numout,9429)
[503]	656	WRITE (numout,9430) ( t2(jpicpt_xad) + t2(jpicpt_yad) + t2(jpicpt_zad) + t2(jpicpt_ldf) + t2(jpicpt_zdf) &
	657	& + t2(jpicpt_npc) + t2(jpicpt_dmp) + t2(jpicpt_qsr) + t2(jpicpt_nsr) ) / tvolt, &
	658	& ( s2(jpicpt_xad) + s2(jpicpt_yad) + s2(jpicpt_zad) + s2(jpicpt_ldf) + s2(jpicpt_zdf) &
	659	& + s2(jpicpt_npc) + s2(jpicpt_dmp) + s2(jpicpt_nsr) ) / tvolt
[215]	660	ENDIF
	661
	662	9420 FORMAT(' tracer**2 trend at it= ', i6, ' : temperature', &
	663	' salinity', /, ' ===============================')
[2528]	664	9421 FORMAT(' zonal advection * t ', e20.13, ' ', e20.13)
	665	9431 FORMAT(' meridional advection * t ', e20.13, ' ', e20.13)
[215]	666	9422 FORMAT(' vertical advection * t ', e20.13, ' ', e20.13)
	667	9423 FORMAT(' horizontal diffusion * t ', e20.13, ' ', e20.13)
	668	9424 FORMAT(' vertical diffusion * t ', e20.13, ' ', e20.13)
[503]	669	9425 FORMAT(' static instability mixing * t ', e20.13, ' ', e20.13)
[215]	670	9426 FORMAT(' damping term * t ', e20.13, ' ', e20.13)
[503]	671	9427 FORMAT(' penetrative qsr * t ', e20.13)
	672	9428 FORMAT(' non solar radiation * t ', e20.13, ' ', e20.13)
[215]	673	9429 FORMAT(' -----------------------------------------------------------------------------')
	674	9430 FORMAT(' total trend t = ', e20.13, ' s = ', e20.13)
	675
	676
	677	IF(lwp) THEN
	678	WRITE (numout,*)
	679	WRITE (numout,*)
	680	WRITE (numout,9440) kt
[503]	681	WRITE (numout,9441) ( tmo(jpicpt_xad)+tmo(jpicpt_yad)+tmo(jpicpt_zad) )/tvolt, &
	682	& ( smo(jpicpt_xad)+smo(jpicpt_yad)+smo(jpicpt_zad) )/tvolt
	683	WRITE (numout,9442) tmo(jpicpt_zl1)/tvolt, smo(jpicpt_zl1)/tvolt
	684	WRITE (numout,9443) tmo(jpicpt_ldf)/tvolt, smo(jpicpt_ldf)/tvolt
	685	WRITE (numout,9444) tmo(jpicpt_zdf)/tvolt, smo(jpicpt_zdf)/tvolt
	686	WRITE (numout,9445) tmo(jpicpt_npc)/tvolt, smo(jpicpt_npc)/tvolt
	687	WRITE (numout,9446) ( t2(jpicpt_xad)+t2(jpicpt_yad)+t2(jpicpt_zad) )/tvolt, &
	688	& ( s2(jpicpt_xad)+s2(jpicpt_yad)+s2(jpicpt_zad) )/tvolt
	689	WRITE (numout,9447) t2(jpicpt_ldf)/tvolt, s2(jpicpt_ldf)/tvolt
	690	WRITE (numout,9448) t2(jpicpt_zdf)/tvolt, s2(jpicpt_zdf)/tvolt
	691	WRITE (numout,9449) t2(jpicpt_npc)/tvolt, s2(jpicpt_npc)/tvolt
[215]	692	ENDIF
	693
	694	9440 FORMAT(' tracer consistency at it= ',i6, &
	695	' : temperature',' salinity',/, &
	696	' ==================================')
[503]	697	9441 FORMAT(' 0 = horizontal+vertical advection + ',e20.13,' ',e20.13)
	698	9442 FORMAT(' 1st lev vertical advection ',e20.13,' ',e20.13)
	699	9443 FORMAT(' 0 = horizontal diffusion ',e20.13,' ',e20.13)
	700	9444 FORMAT(' 0 = vertical diffusion ',e20.13,' ',e20.13)
	701	9445 FORMAT(' 0 = static instability mixing ',e20.13,' ',e20.13)
	702	9446 FORMAT(' 0 = horizontal+vertical advection * t ',e20.13,' ',e20.13)
	703	9447 FORMAT(' 0 > horizontal diffusion * t ',e20.13,' ',e20.13)
	704	9448 FORMAT(' 0 > vertical diffusion * t ',e20.13,' ',e20.13)
	705	9449 FORMAT(' 0 > static instability mixing * t ',e20.13,' ',e20.13)
	706	!
[215]	707	ENDIF
[503]	708	!
[215]	709	END SUBROUTINE trd_twr
	710
	711	# else
	712	!!----------------------------------------------------------------------
	713	!! Default case : Empty module
	714	!!----------------------------------------------------------------------
[601]	715	INTERFACE trd_icp
	716	MODULE PROCEDURE trd_2d, trd_3d
	717	END INTERFACE
	718
[215]	719	CONTAINS
[2528]	720	SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype ) ! Empty routine
[503]	721	REAL, DIMENSION(:,:) :: ptrd2dx, ptrd2dy
[601]	722	INTEGER , INTENT(in ) :: ktrd ! tracer trend index
	723	CHARACTER(len=3) , INTENT(in ) :: ctype ! momentum ('DYN') or tracers ('TRA') trends
[521]	724	WRITE(,) 'trd_2d: You should not have seen this print! error ?', &
[2528]	725	& ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype
[215]	726	END SUBROUTINE trd_2d
[2528]	727	SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype ) ! Empty routine
[503]	728	REAL, DIMENSION(:,:,:) :: ptrd3dx, ptrd3dy
[601]	729	INTEGER , INTENT(in ) :: ktrd ! tracer trend index
	730	CHARACTER(len=3) , INTENT(in ) :: ctype ! momentum ('DYN') or tracers ('TRA') trends
[521]	731	WRITE(,) 'trd_3d: You should not have seen this print! error ?', &
[2528]	732	& ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype
[215]	733	END SUBROUTINE trd_3d
	734	SUBROUTINE trd_icp_init ! Empty routine
	735	END SUBROUTINE trd_icp_init
	736	SUBROUTINE trd_dwr( kt ) ! Empty routine
	737	WRITE(,) 'trd_dwr: You should not have seen this print! error ?', kt
	738	END SUBROUTINE trd_dwr
	739	SUBROUTINE trd_twr( kt ) ! Empty routine
	740	WRITE(,) 'trd_twr: You should not have seen this print! error ?', kt
	741	END SUBROUTINE trd_twr
	742	#endif
	743
	744	!!======================================================================
	745	END MODULE trdicp

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

Download in other formats: