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/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

Context Navigation

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90 @ 3316

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

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

Download in other formats: