Context Navigation

diawri.F90

Last change on this file was 15260, checked in by acc, 3 years ago
Untested, initial port of GEOMETRIC changes. Includes a bug fix to eken diagnostic in diawri.F90 which, if confirmed, is also relevant to current 4.0-HEAD. #2722
Property svn:keywords set to `Id`
File size: 52.1 KB

Rev	Line
[9403]	1	MODULE diawri
	2	!!======================================================================
	3	!! * MODULE diawri *
	4	!! Ocean diagnostics : write ocean output files
	5	!!=====================================================================
	6	!! History : OPA ! 1991-03 (M.-A. Foujols) Original code
	7	!! 4.0 ! 1991-11 (G. Madec)
	8	!! ! 1992-06 (M. Imbard) correction restart file
	9	!! ! 1992-07 (M. Imbard) split into diawri and rstwri
	10	!! ! 1993-03 (M. Imbard) suppress writibm
	11	!! ! 1998-01 (C. Levy) NETCDF format using ioipsl INTERFACE
	12	!! ! 1999-02 (E. Guilyardi) name of netCDF files + variables
	13	!! 8.2 ! 2000-06 (M. Imbard) Original code (diabort.F)
	14	!! NEMO 1.0 ! 2002-06 (A.Bozec, E. Durand) Original code (diainit.F)
	15	!! - ! 2002-09 (G. Madec) F90: Free form and module
	16	!! - ! 2002-12 (G. Madec) merge of diabort and diainit, F90
	17	!! ! 2005-11 (V. Garnier) Surface pressure gradient organization
	18	!! 3.2 ! 2008-11 (B. Lemaire) creation from old diawri
	19	!! 3.7 ! 2014-01 (G. Madec) remove eddy induced velocity from no-IOM output
	20	!! ! change name of output variables in dia_wri_state
	21	!!----------------------------------------------------------------------
	22
	23	!!----------------------------------------------------------------------
	24	!! dia_wri : create the standart output files
	25	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
	26	!!----------------------------------------------------------------------
	27	USE oce ! ocean dynamics and tracers
	28	USE dom_oce ! ocean space and time domain
	29	USE phycst ! physical constants
	30	USE dianam ! build name of file (routine)
	31	USE diahth ! thermocline diagnostics
	32	USE dynadv , ONLY: ln_dynadv_vec
	33	USE icb_oce ! Icebergs
	34	USE icbdia ! Iceberg budgets
	35	USE ldftra ! lateral physics: eddy diffusivity coef.
	36	USE ldfdyn ! lateral physics: eddy viscosity coef.
	37	USE sbc_oce ! Surface boundary condition: ocean fields
	38	USE sbc_ice ! Surface boundary condition: ice fields
	39	USE sbcssr ! restoring term toward SST/SSS climatology
	40	USE sbcwave ! wave parameters
	41	USE wet_dry ! wetting and drying
	42	USE zdf_oce ! ocean vertical physics
	43	USE zdfdrg ! ocean vertical physics: top/bottom friction
	44	USE zdfmxl ! mixed layer
	45	!
	46	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	47	USE in_out_manager ! I/O manager
	48	USE dia25h ! 25h Mean output
	49	USE iom !
	50	USE ioipsl !
	51
[9576]	52	#if defined key_si3
[10425]	53	USE ice
[9403]	54	USE icewri
	55	#endif
	56	USE lib_mpp ! MPP library
	57	USE timing ! preformance summary
	58	USE diurnal_bulk ! diurnal warm layer
	59	USE cool_skin ! Cool skin
	60
	61	IMPLICIT NONE
	62	PRIVATE
	63
	64	PUBLIC dia_wri ! routines called by step.F90
	65	PUBLIC dia_wri_state
	66	PUBLIC dia_wri_alloc ! Called by nemogcm module
	67
	68	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
	69	INTEGER :: nb_T , ndim_bT ! grid_T file
	70	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	71	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	72	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
	73	INTEGER :: ndex(1) ! ???
	74	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
	75	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
	76	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
	77
	78	!! * Substitutions
	79	# include "vectopt_loop_substitute.h90"
	80	!!----------------------------------------------------------------------
[10073]	81	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10074]	82	!! $Id$
[10073]	83	!! Software governed by the CeCILL license (see ./LICENSE)
[9403]	84	!!----------------------------------------------------------------------
	85	CONTAINS
	86
	87	#if defined key_iomput
	88	!!----------------------------------------------------------------------
	89	!! 'key_iomput' use IOM library
	90	!!----------------------------------------------------------------------
[9652]	91	INTEGER FUNCTION dia_wri_alloc()
	92	!
	93	dia_wri_alloc = 0
	94	!
	95	END FUNCTION dia_wri_alloc
	96
[10425]	97
[9403]	98	SUBROUTINE dia_wri( kt )
	99	!!---------------------------------------------------------------------
	100	!! * ROUTINE dia_wri *
	101	!!
	102	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	103	!! NETCDF format is used by default
	104	!!
	105	!! ** Method : use iom_put
	106	!!----------------------------------------------------------------------
	107	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	108	!!
	109	INTEGER :: ji, jj, jk ! dummy loop indices
	110	INTEGER :: ikbot ! local integer
	111	REAL(wp):: zztmp , zztmpx ! local scalar
	112	REAL(wp):: zztmp2, zztmpy ! - -
	113	REAL(wp):: ze3 ! - -
	114	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
	115	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
	116	REAL(wp), DIMENSION(jpi,jpj,jpk) :: bu, bv ! volume of u- and v-boxes
	117	REAL(wp), DIMENSION(jpi,jpj,jpk) :: r1_bt ! inverse of t-box volume
	118	!!----------------------------------------------------------------------
	119	!
	120	IF( ln_timing ) CALL timing_start('dia_wri')
	121	!
	122	! Output the initial state and forcings
	123	IF( ninist == 1 ) THEN
[10425]	124	CALL dia_wri_state( 'output.init' )
[9403]	125	ninist = 0
	126	ENDIF
	127
	128	! Output of initial vertical scale factor
	129	CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10425]	130	CALL iom_put("e3u_0", e3u_0(:,:,:) )
	131	CALL iom_put("e3v_0", e3v_0(:,:,:) )
[9403]	132	!
	133	CALL iom_put( "e3t" , e3t_n(:,:,:) )
	134	CALL iom_put( "e3u" , e3u_n(:,:,:) )
	135	CALL iom_put( "e3v" , e3v_n(:,:,:) )
	136	CALL iom_put( "e3w" , e3w_n(:,:,:) )
	137	IF( iom_use("e3tdef") ) &
	138	CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
	139
	140	IF( ll_wd ) THEN
	141	CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)
	142	ELSE
	143	CALL iom_put( "ssh" , sshn ) ! sea surface height
	144	ENDIF
	145
	146	IF( iom_use("wetdep") ) & ! wet depth
	147	CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
	148
	149	CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature
	150	CALL iom_put( "sst", tsn(:,:,1,jp_tem) ) ! surface temperature
	151	IF ( iom_use("sbt") ) THEN
	152	DO jj = 1, jpj
	153	DO ji = 1, jpi
	154	ikbot = mbkt(ji,jj)
	155	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
	156	END DO
	157	END DO
	158	CALL iom_put( "sbt", z2d ) ! bottom temperature
	159	ENDIF
	160
	161	CALL iom_put( "soce", tsn(:,:,:,jp_sal) ) ! 3D salinity
	162	CALL iom_put( "sss", tsn(:,:,1,jp_sal) ) ! surface salinity
	163	IF ( iom_use("sbs") ) THEN
	164	DO jj = 1, jpj
	165	DO ji = 1, jpi
	166	ikbot = mbkt(ji,jj)
	167	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
	168	END DO
	169	END DO
	170	CALL iom_put( "sbs", z2d ) ! bottom salinity
	171	ENDIF
	172
	173	IF ( iom_use("taubot") ) THEN ! bottom stress
	174	zztmp = rau0 * 0.25
	175	z2d(:,:) = 0._wp
	176	DO jj = 2, jpjm1
	177	DO ji = fs_2, fs_jpim1 ! vector opt.
	178	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * un(ji ,jj,mbku(ji ,jj)) )**2 &
	179	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj)) )**2 &
	180	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vn(ji,jj ,mbkv(ji,jj )) )**2 &
	181	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1)) )**2
	182	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
	183	!
	184	END DO
	185	END DO
[10425]	186	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[9403]	187	CALL iom_put( "taubot", z2d )
	188	ENDIF
	189
	190	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
	191	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
	192	IF ( iom_use("sbu") ) THEN
	193	DO jj = 1, jpj
	194	DO ji = 1, jpi
	195	ikbot = mbku(ji,jj)
	196	z2d(ji,jj) = un(ji,jj,ikbot)
	197	END DO
	198	END DO
	199	CALL iom_put( "sbu", z2d ) ! bottom i-current
	200	ENDIF
	201
	202	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
	203	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
	204	IF ( iom_use("sbv") ) THEN
	205	DO jj = 1, jpj
	206	DO ji = 1, jpi
	207	ikbot = mbkv(ji,jj)
	208	z2d(ji,jj) = vn(ji,jj,ikbot)
	209	END DO
	210	END DO
	211	CALL iom_put( "sbv", z2d ) ! bottom j-current
	212	ENDIF
	213
[15260]	214	! ! vertical velocity
	215	IF( ln_zad_Aimp ) THEN ; CALL iom_put( "woce", wn + wi ) ! explicit plus implicit parts
	216	ELSE ; CALL iom_put( "woce", wn )
	217	ENDIF
	218
[9403]	219	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
[15260]	220	! ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[9403]	221	DO jk = 1, jpk
[15260]	222	IF( ln_zad_Aimp ) THEN
	223	z3d(:,:,jk) = rau0 * e1e2t(:,:) * ( wn(:,:,jk) + wi(:,:,jk) )
	224	ELSE
	225	z3d(:,:,jk) = rau0 * e1e2t(:,:) * wn(:,:,jk)
	226	ENDIF
[9403]	227	END DO
	228	CALL iom_put( "w_masstr" , z3d )
[15260]	229	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d * z3d )
[9403]	230	ENDIF
	231
	232	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	233	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
	234	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
	235
	236	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
	237	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
	238
[13284]	239	IF ( iom_use("socegrad") .OR. iom_use("socegrad2") ) THEN
[9403]	240	z3d(:,:,jpk) = 0.
	241	DO jk = 1, jpkm1
	242	DO jj = 2, jpjm1 ! sal gradient
	243	DO ji = fs_2, fs_jpim1 ! vector opt.
	244	zztmp = tsn(ji,jj,jk,jp_sal)
	245	zztmpx = ( tsn(ji+1,jj,jk,jp_sal) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,jk,jp_sal) ) * r1_e1u(ji-1,jj)
	246	zztmpy = ( tsn(ji,jj+1,jk,jp_sal) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,jk,jp_sal) ) * r1_e2v(ji,jj-1)
	247	z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	248	& * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk)
	249	END DO
	250	END DO
	251	END DO
[10425]	252	CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
[13284]	253	CALL iom_put( "socegrad2", z3d ) ! square of module of sal gradient
[9403]	254	z3d(:,:,:) = SQRT( z3d(:,:,:) )
[13284]	255	CALL iom_put( "socegrad" , z3d ) ! module of sal gradient
[9403]	256	ENDIF
	257
	258	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
	259	DO jj = 2, jpjm1 ! sst gradient
	260	DO ji = fs_2, fs_jpim1 ! vector opt.
	261	zztmp = tsn(ji,jj,1,jp_tem)
	262	zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,1,jp_tem) ) * r1_e1u(ji-1,jj)
	263	zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
	264	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	265	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
	266	END DO
	267	END DO
[10425]	268	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[9403]	269	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
	270	z2d(:,:) = SQRT( z2d(:,:) )
	271	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
	272	ENDIF
	273
	274	! heat and salt contents
	275	IF( iom_use("heatc") ) THEN
	276	z2d(:,:) = 0._wp
	277	DO jk = 1, jpkm1
	278	DO jj = 1, jpj
	279	DO ji = 1, jpi
	280	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
	281	END DO
	282	END DO
	283	END DO
	284	CALL iom_put( "heatc", rau0_rcp * z2d ) ! vertically integrated heat content (J/m2)
	285	ENDIF
	286
	287	IF( iom_use("saltc") ) THEN
	288	z2d(:,:) = 0._wp
	289	DO jk = 1, jpkm1
	290	DO jj = 1, jpj
	291	DO ji = 1, jpi
	292	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
	293	END DO
	294	END DO
	295	END DO
	296	CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
	297	ENDIF
	298	!
	299	IF( iom_use("salt2c") ) THEN
	300	z2d(:,:) = 0._wp
	301	DO jk = 1, jpkm1
	302	DO jj = 1, jpj
	303	DO ji = 1, jpi
	304	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
	305	END DO
	306	END DO
	307	END DO
[13284]	308	CALL iom_put( "salt2c", rau0 * z2d ) ! vertically integrated squared salt content (PSU*kg/m2)
[9403]	309	ENDIF
	310	!
[13284]	311	IF ( iom_use("eken") .OR. iom_use("eken_int") ) THEN
[9403]	312	z3d(:,:,jpk) = 0._wp
	313	DO jk = 1, jpkm1
[13284]	314	DO jj = 2, jpjm1
	315	DO ji = 2, jpim1
[9403]	316	zztmpx = 0.5 * ( un(ji-1,jj ,jk) + un(ji,jj,jk) )
	317	zztmpy = 0.5 * ( vn(ji ,jj-1,jk) + vn(ji,jj,jk) )
	318	z3d(ji,jj,jk) = 0.5 * ( zztmpxzztmpx + zztmpyzztmpy )
	319	END DO
	320	END DO
	321	END DO
[10425]	322	CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
[9403]	323	CALL iom_put( "eken", z3d ) ! kinetic energy
	324
	325	z2d(:,:) = 0._wp
	326	DO jk = 1, jpkm1
	327	DO jj = 1, jpj
	328	DO ji = 1, jpi
[13284]	329	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
[9403]	330	END DO
	331	END DO
	332	END DO
[13284]	333	CALL iom_put( "eken_int", z2d ) ! vertically integrated kinetic energy
[9403]	334	ENDIF
	335	!
	336	CALL iom_put( "hdiv", hdivn ) ! Horizontal divergence
	337
	338	IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN
	339
	340	z3d(:,:,jpk) = 0._wp
	341	DO jk = 1, jpkm1
	342	DO jj = 1, jpjm1
	343	DO ji = 1, fs_jpim1 ! vector opt.
	344	z3d(ji,jj,jk) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &
[13284]	345	& - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj)
[9403]	346	END DO
	347	END DO
	348	END DO
[10425]	349	CALL lbc_lnk( 'diawri', z3d, 'F', 1. )
[9403]	350	CALL iom_put( "relvor", z3d ) ! relative vorticity
	351
	352	DO jk = 1, jpkm1
	353	DO jj = 1, jpj
	354	DO ji = 1, jpi
	355	z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)
	356	END DO
	357	END DO
	358	END DO
	359	CALL iom_put( "absvor", z3d ) ! absolute vorticity
	360
	361	DO jk = 1, jpkm1
	362	DO jj = 1, jpjm1
	363	DO ji = 1, fs_jpim1 ! vector opt.
	364	ze3 = ( e3t_n(ji,jj+1,jk)tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
	365	& + e3t_n(ji,jj ,jk)tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)tmask(ji+1,jj ,jk) )
	366	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	367	ELSE ; ze3 = 0._wp
	368	ENDIF
	369	z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)
	370	END DO
	371	END DO
	372	END DO
[10425]	373	CALL lbc_lnk( 'diawri', z3d, 'F', 1. )
[9403]	374	CALL iom_put( "potvor", z3d ) ! potential vorticity
	375
	376	ENDIF
	377
	378	!
	379	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
	380	z3d(:,:,jpk) = 0.e0
	381	z2d(:,:) = 0.e0
	382	DO jk = 1, jpkm1
	383	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
	384	z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
	385	END DO
	386	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
	387	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
	388	ENDIF
	389
	390	IF( iom_use("u_heattr") ) THEN
	391	z2d(:,:) = 0._wp
	392	DO jk = 1, jpkm1
	393	DO jj = 2, jpjm1
	394	DO ji = fs_2, fs_jpim1 ! vector opt.
	395	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
	396	END DO
	397	END DO
	398	END DO
[10425]	399	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9403]	400	CALL iom_put( "u_heattr", 0.5rcp z2d ) ! heat transport in i-direction
	401	ENDIF
	402
	403	IF( iom_use("u_salttr") ) THEN
	404	z2d(:,:) = 0.e0
	405	DO jk = 1, jpkm1
	406	DO jj = 2, jpjm1
	407	DO ji = fs_2, fs_jpim1 ! vector opt.
	408	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
	409	END DO
	410	END DO
	411	END DO
[10425]	412	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9403]	413	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
	414	ENDIF
	415
	416
	417	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
	418	z3d(:,:,jpk) = 0.e0
	419	DO jk = 1, jpkm1
	420	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
	421	END DO
	422	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
	423	ENDIF
	424
	425	IF( iom_use("v_heattr") ) THEN
	426	z2d(:,:) = 0.e0
	427	DO jk = 1, jpkm1
	428	DO jj = 2, jpjm1
	429	DO ji = fs_2, fs_jpim1 ! vector opt.
	430	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
	431	END DO
	432	END DO
	433	END DO
[10425]	434	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9403]	435	CALL iom_put( "v_heattr", 0.5rcp z2d ) ! heat transport in j-direction
	436	ENDIF
	437
	438	IF( iom_use("v_salttr") ) THEN
	439	z2d(:,:) = 0._wp
	440	DO jk = 1, jpkm1
	441	DO jj = 2, jpjm1
	442	DO ji = fs_2, fs_jpim1 ! vector opt.
	443	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
	444	END DO
	445	END DO
	446	END DO
[10425]	447	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9403]	448	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
	449	ENDIF
	450
	451	IF( iom_use("tosmint") ) THEN
	452	z2d(:,:) = 0._wp
	453	DO jk = 1, jpkm1
	454	DO jj = 2, jpjm1
	455	DO ji = fs_2, fs_jpim1 ! vector opt.
	456	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
	457	END DO
	458	END DO
	459	END DO
[10425]	460	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9403]	461	CALL iom_put( "tosmint", rau0 * z2d ) ! Vertical integral of temperature
	462	ENDIF
	463	IF( iom_use("somint") ) THEN
	464	z2d(:,:)=0._wp
	465	DO jk = 1, jpkm1
	466	DO jj = 2, jpjm1
	467	DO ji = fs_2, fs_jpim1 ! vector opt.
	468	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
	469	END DO
	470	END DO
	471	END DO
[10425]	472	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9403]	473	CALL iom_put( "somint", rau0 * z2d ) ! Vertical integral of salinity
	474	ENDIF
	475
	476	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
	477	!
	478
	479	IF (ln_dia25h) CALL dia_25h( kt ) ! 25h averaging
	480
	481	IF( ln_timing ) CALL timing_stop('dia_wri')
	482	!
	483	END SUBROUTINE dia_wri
	484
	485	#else
	486	!!----------------------------------------------------------------------
	487	!! Default option use IOIPSL library
	488	!!----------------------------------------------------------------------
[10425]	489
[9652]	490	INTEGER FUNCTION dia_wri_alloc()
	491	!!----------------------------------------------------------------------
	492	INTEGER, DIMENSION(2) :: ierr
	493	!!----------------------------------------------------------------------
	494	ierr = 0
	495	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	496	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	497	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
	498	!
	499	dia_wri_alloc = MAXVAL(ierr)
[10425]	500	CALL mpp_sum( 'diawri', dia_wri_alloc )
[9652]	501	!
	502	END FUNCTION dia_wri_alloc
[10425]	503
	504
[9403]	505	SUBROUTINE dia_wri( kt )
	506	!!---------------------------------------------------------------------
	507	!! * ROUTINE dia_wri *
	508	!!
	509	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	510	!! NETCDF format is used by default
	511	!!
	512	!! ** Method : At the beginning of the first time step (nit000),
	513	!! define all the NETCDF files and fields
	514	!! At each time step call histdef to compute the mean if ncessary
[11536]	515	!! Each nn_write time step, output the instantaneous or mean fields
[9403]	516	!!----------------------------------------------------------------------
	517	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	518	!
	519	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	520	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	521	INTEGER :: inum = 11 ! temporary logical unit
	522	INTEGER :: ji, jj, jk ! dummy loop indices
	523	INTEGER :: ierr ! error code return from allocation
	524	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
	525	INTEGER :: jn, ierror ! local integers
	526	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
	527	!
	528	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
	529	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
	530	!!----------------------------------------------------------------------
	531	!
	532	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[10425]	533	CALL dia_wri_state( 'output.init' )
[9403]	534	ninist = 0
	535	ENDIF
	536	!
[11536]	537	IF( nn_write == -1 ) RETURN ! we will never do any output
	538	!
	539	IF( ln_timing ) CALL timing_start('dia_wri')
	540	!
[9403]	541	! 0. Initialisation
	542	! -----------------
	543
	544	ll_print = .FALSE. ! local variable for debugging
	545	ll_print = ll_print .AND. lwp
	546
	547	! Define frequency of output and means
	548	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
	549	#if defined key_diainstant
[11536]	550	zsto = nn_write * rdt
[9403]	551	clop = "inst("//TRIM(clop)//")"
	552	#else
	553	zsto=rdt
	554	clop = "ave("//TRIM(clop)//")"
	555	#endif
[11536]	556	zout = nn_write * rdt
[9403]	557	zmax = ( nitend - nit000 + 1 ) * rdt
	558
	559	! Define indices of the horizontal output zoom and vertical limit storage
	560	iimi = 1 ; iima = jpi
	561	ijmi = 1 ; ijma = jpj
	562	ipk = jpk
	563
	564	! define time axis
	565	it = kt
	566	itmod = kt - nit000 + 1
	567
	568
	569	! 1. Define NETCDF files and fields at beginning of first time step
	570	! -----------------------------------------------------------------
	571
	572	IF( kt == nit000 ) THEN
	573
	574	! Define the NETCDF files (one per grid)
	575
	576	! Compute julian date from starting date of the run
	577	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	578	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
	579	IF(lwp)WRITE(numout,*)
	580	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	581	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	582	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	583	' limit storage in depth = ', ipk
	584
	585	! WRITE root name in date.file for use by postpro
	586	IF(lwp) THEN
[11536]	587	CALL dia_nam( clhstnam, nn_write,' ' )
[9403]	588	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
	589	WRITE(inum,*) clhstnam
	590	CLOSE(inum)
	591	ENDIF
	592
	593	! Define the T grid FILE ( nid_T )
	594
[11536]	595	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[9403]	596	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	597	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	598	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	599	& nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
	600	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
	601	& "m", ipk, gdept_1d, nz_T, "down" )
	602	! ! Index of ocean points
	603	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	604	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
	605	!
	606	IF( ln_icebergs ) THEN
	607	!
	608	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	609	!! that routine is called from nemogcm, so do it here immediately before its needed
	610	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	611	CALL mpp_sum( 'diawri', ierror )
[9403]	612	IF( ierror /= 0 ) THEN
	613	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	614	RETURN
	615	ENDIF
	616	!
	617	!! iceberg vertical coordinate is class number
	618	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	619	& "number", nclasses, class_num, nb_T )
	620	!
	621	!! each class just needs the surface index pattern
	622	ndim_bT = 3
	623	DO jn = 1,nclasses
	624	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	625	ENDDO
	626	!
	627	ENDIF
	628
	629	! Define the U grid FILE ( nid_U )
	630
[11536]	631	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[9403]	632	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	633	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	634	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	635	& nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
	636	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
	637	& "m", ipk, gdept_1d, nz_U, "down" )
	638	! ! Index of ocean points
	639	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	640	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	641
	642	! Define the V grid FILE ( nid_V )
	643
[11536]	644	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[9403]	645	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	646	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	647	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	648	& nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
	649	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
	650	& "m", ipk, gdept_1d, nz_V, "down" )
	651	! ! Index of ocean points
	652	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	653	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	654
	655	! Define the W grid FILE ( nid_W )
	656
[11536]	657	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[9403]	658	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	659	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	660	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	661	& nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
	662	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
	663	& "m", ipk, gdepw_1d, nz_W, "down" )
	664
	665
	666	! Declare all the output fields as NETCDF variables
	667
	668	! !!! nid_T : 3D
	669	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	670	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	671	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	672	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	673	IF( .NOT.ln_linssh ) THEN
	674	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	675	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	676	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	677	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	678	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	679	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	680	ENDIF
	681	! !!! nid_T : 2D
	682	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	683	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	684	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	685	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	686	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
	687	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	688	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
	689	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	690	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	691	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	692	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
	693	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	694	IF( ln_linssh ) THEN
	695	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
	696	& , "KgC/m2/s", & ! sosst_cd
	697	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	698	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
	699	& , "KgPSU/m2/s",& ! sosss_cd
	700	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	701	ENDIF
	702	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
	703	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	704	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	705	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	706	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	707	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	708	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	709	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	710	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
	711	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	712	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	713	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	714	!
	715	IF( ln_icebergs ) THEN
	716	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	717	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	718	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	719	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	720	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	721	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	722	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	723	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	724	IF( ln_bergdia ) THEN
	725	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	726	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	727	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	728	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	729	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	730	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	731	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	732	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	733	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	734	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	735	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	736	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	737	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	738	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	739	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	740	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	741	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	742	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	743	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	744	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	745	ENDIF
	746	ENDIF
	747
[11536]	748	IF( ln_ssr ) THEN
[9403]	749	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	750	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	751	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	752	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	753	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	754	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	755	ENDIF
[11536]	756
[9403]	757	clmx ="l_max(only(x))" ! max index on a period
	758	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	759	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
	760	#if defined key_diahth
	761	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	762	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	763	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	764	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	765	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	766	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	767	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
	768	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	769	#endif
	770
	771	CALL histend( nid_T, snc4chunks=snc4set )
	772
	773	! !!! nid_U : 3D
	774	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	775	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	776	IF( ln_wave .AND. ln_sdw) THEN
	777	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	778	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	779	ENDIF
	780	! !!! nid_U : 2D
	781	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
	782	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	783
	784	CALL histend( nid_U, snc4chunks=snc4set )
	785
	786	! !!! nid_V : 3D
	787	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	788	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	789	IF( ln_wave .AND. ln_sdw) THEN
	790	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	791	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	792	ENDIF
	793	! !!! nid_V : 2D
	794	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
	795	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	796
	797	CALL histend( nid_V, snc4chunks=snc4set )
	798
	799	! !!! nid_W : 3D
	800	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	801	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	802	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	803	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	804	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
	805	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	806
	807	IF( ln_zdfddm ) THEN
	808	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	809	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	810	ENDIF
	811
	812	IF( ln_wave .AND. ln_sdw) THEN
	813	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	814	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	815	ENDIF
	816	! !!! nid_W : 2D
	817	CALL histend( nid_W, snc4chunks=snc4set )
	818
	819	IF(lwp) WRITE(numout,*)
	820	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	821	IF(ll_print) CALL FLUSH(numout )
	822
	823	ENDIF
	824
	825	! 2. Start writing data
	826	! ---------------------
	827
	828	! ndex(1) est utilise ssi l'avant dernier argument est different de
	829	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	830	! donne le nombre d'elements, et ndex la liste des indices a sortir
	831
[11536]	832	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
[9403]	833	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	834	WRITE(numout,*) '~~~~~~ '
	835	ENDIF
	836
	837	IF( .NOT.ln_linssh ) THEN
	838	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	839	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	840	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	841	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
	842	ELSE
	843	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	844	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	845	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	846	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	847	ENDIF
	848	IF( .NOT.ln_linssh ) THEN
	849	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
	850	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	851	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
	852	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	853	ENDIF
	854	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
	855	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
	856	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
	857	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	858	! (includes virtual salt flux beneath ice
	859	! in linear free surface case)
	860	IF( ln_linssh ) THEN
	861	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
	862	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
	863	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
	864	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	865	ENDIF
	866	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
	867	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
	868	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
	869	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
	870	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
	871	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
	872	!
	873	IF( ln_icebergs ) THEN
	874	!
	875	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	876	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	877	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	878	!
	879	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	880	!
	881	IF( ln_bergdia ) THEN
	882	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	883	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	884	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	885	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	886	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	887	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	888	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	889	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	890	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	891	!
	892	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	893	ENDIF
	894	ENDIF
	895
[11536]	896	IF( ln_ssr ) THEN
[9403]	897	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	898	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[11536]	899	zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[9403]	900	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	901	ENDIF
	902	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	903	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
	904
	905	#if defined key_diahth
	906	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	907	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	908	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	909	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	910	#endif
	911
	912	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
	913	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
	914
	915	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
	916	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
	917
	918	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
	919	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
	920	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	921	IF( ln_zdfddm ) THEN
	922	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
	923	ENDIF
	924
	925	IF( ln_wave .AND. ln_sdw ) THEN
	926	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	927	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	928	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
	929	ENDIF
	930
	931	! 3. Close all files
	932	! ---------------------------------------
	933	IF( kt == nitend ) THEN
	934	CALL histclo( nid_T )
	935	CALL histclo( nid_U )
	936	CALL histclo( nid_V )
	937	CALL histclo( nid_W )
	938	ENDIF
	939	!
	940	IF( ln_timing ) CALL timing_stop('dia_wri')
	941	!
	942	END SUBROUTINE dia_wri
	943	#endif
	944
[10425]	945	SUBROUTINE dia_wri_state( cdfile_name )
[9403]	946	!!---------------------------------------------------------------------
	947	!! * ROUTINE dia_wri_state *
	948	!!
	949	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	950	!! the instantaneous ocean state and forcing fields.
	951	!! Used to find errors in the initial state or save the last
	952	!! ocean state in case of abnormal end of a simulation
	953	!!
	954	!! ** Method : NetCDF files using ioipsl
	955	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	956	!! File 'output.abort.nc' is created in case of abnormal job end
	957	!!----------------------------------------------------------------------
	958	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	959	!!
	960	INTEGER :: inum
[9403]	961	!!----------------------------------------------------------------------
	962	!
	963	IF(lwp) WRITE(numout,*)
	964	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	965	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
[10425]	966	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
[9403]	967
[9576]	968	#if defined key_si3
[10425]	969	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
[9403]	970	#else
[10425]	971	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
[9403]	972	#endif
	973
[10425]	974	CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) ) ! now temperature
	975	CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) ) ! now salinity
	976	CALL iom_rstput( 0, 0, inum, 'sossheig', sshn ) ! sea surface height
	977	CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity
	978	CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity
	979	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity
[9403]	980	IF( ALLOCATED(ahtu) ) THEN
[10425]	981	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	982	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[9403]	983	ENDIF
	984	IF( ALLOCATED(ahmt) ) THEN
[10425]	985	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	986	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[9403]	987	ENDIF
[10425]	988	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	989	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	990	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	991	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	992	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	993	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[9403]	994	IF( .NOT.ln_linssh ) THEN
[10425]	995	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
	996	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
	997	END IF
[9403]	998	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	999	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	1000	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	1001	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[9403]	1002	ENDIF
[10425]	1003
	1004	#if defined key_si3
	1005	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
	1006	CALL ice_wri_state( inum )
[9403]	1007	ENDIF
	1008	#endif
[10425]	1009	!
	1010	CALL iom_close( inum )
[9403]	1011	!
	1012	END SUBROUTINE dia_wri_state
	1013
	1014	!!======================================================================
	1015	END MODULE diawri

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

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

Context Navigation

source: NEMO/branches/NERC/dev_r4.0.6_GEOMETRIC/tests/CANAL/MY_SRC/diawri.F90

Download in other formats: