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.

diawri.F90 in NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/tests/CANAL/MY_SRC – NEMO

Context Navigation

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/tests/CANAL/MY_SRC/diawri.F90 @ 12756

Last change on this file since 12756 was 12756, checked in by clem, 4 years ago
fix xml again
Property svn:keywords set to `Id`
File size: 51.8 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
	214	CALL iom_put( "woce", wn ) ! vertical velocity
	215	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
	216	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
	217	z2d(:,:) = rau0 * e1e2t(:,:)
	218	DO jk = 1, jpk
	219	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
	220	END DO
	221	CALL iom_put( "w_masstr" , z3d )
	222	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
	223	ENDIF
	224
	225	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	226	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
	227	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
	228
	229	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
	230	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
	231
[12756]	232	IF ( iom_use("socegrad") .OR. iom_use("socegrad2") ) THEN
[9403]	233	z3d(:,:,jpk) = 0.
	234	DO jk = 1, jpkm1
	235	DO jj = 2, jpjm1 ! sal gradient
	236	DO ji = fs_2, fs_jpim1 ! vector opt.
	237	zztmp = tsn(ji,jj,jk,jp_sal)
	238	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)
	239	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)
	240	z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	241	& * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk)
	242	END DO
	243	END DO
	244	END DO
[10425]	245	CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
[12756]	246	CALL iom_put( "socegrad2", z3d ) ! square of module of sal gradient
[9403]	247	z3d(:,:,:) = SQRT( z3d(:,:,:) )
[12756]	248	CALL iom_put( "socegrad" , z3d ) ! module of sal gradient
[9403]	249	ENDIF
	250
	251	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
	252	DO jj = 2, jpjm1 ! sst gradient
	253	DO ji = fs_2, fs_jpim1 ! vector opt.
	254	zztmp = tsn(ji,jj,1,jp_tem)
	255	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)
	256	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)
	257	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	258	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
	259	END DO
	260	END DO
[10425]	261	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[9403]	262	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
	263	z2d(:,:) = SQRT( z2d(:,:) )
	264	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
	265	ENDIF
	266
	267	! heat and salt contents
	268	IF( iom_use("heatc") ) THEN
	269	z2d(:,:) = 0._wp
	270	DO jk = 1, jpkm1
	271	DO jj = 1, jpj
	272	DO ji = 1, jpi
	273	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
	274	END DO
	275	END DO
	276	END DO
	277	CALL iom_put( "heatc", rau0_rcp * z2d ) ! vertically integrated heat content (J/m2)
	278	ENDIF
	279
	280	IF( iom_use("saltc") ) THEN
	281	z2d(:,:) = 0._wp
	282	DO jk = 1, jpkm1
	283	DO jj = 1, jpj
	284	DO ji = 1, jpi
	285	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
	286	END DO
	287	END DO
	288	END DO
	289	CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
	290	ENDIF
	291	!
	292	IF( iom_use("salt2c") ) THEN
	293	z2d(:,:) = 0._wp
	294	DO jk = 1, jpkm1
	295	DO jj = 1, jpj
	296	DO ji = 1, jpi
	297	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)
	298	END DO
	299	END DO
	300	END DO
[12756]	301	CALL iom_put( "salt2c", rau0 * z2d ) ! vertically integrated squared salt content (PSU*kg/m2)
[9403]	302	ENDIF
	303	!
[12756]	304	IF ( iom_use("eken") .OR. iom_use("eken_int") ) THEN
[9403]	305	z3d(:,:,jpk) = 0._wp
	306	DO jk = 1, jpkm1
[12756]	307	DO jj = 2, jpjm1
	308	DO ji = 2, jpim1
[9403]	309	zztmpx = 0.5 * ( un(ji-1,jj ,jk) + un(ji,jj,jk) )
	310	zztmpy = 0.5 * ( vn(ji ,jj-1,jk) + vn(ji,jj,jk) )
	311	z3d(ji,jj,jk) = 0.5 * ( zztmpxzztmpx + zztmpyzztmpy )
	312	END DO
	313	END DO
	314	END DO
[10425]	315	CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
[9403]	316	CALL iom_put( "eken", z3d ) ! kinetic energy
	317
	318	z2d(:,:) = 0._wp
	319	DO jk = 1, jpkm1
	320	DO jj = 1, jpj
	321	DO ji = 1, jpi
[12756]	322	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
[9403]	323	END DO
	324	END DO
	325	END DO
[12756]	326	CALL iom_put( "eken_int", z2d ) ! vertically integrated kinetic energy
[9403]	327	ENDIF
	328	!
	329	CALL iom_put( "hdiv", hdivn ) ! Horizontal divergence
	330
	331	IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN
	332
	333	z3d(:,:,jpk) = 0._wp
	334	DO jk = 1, jpkm1
	335	DO jj = 1, jpjm1
	336	DO ji = 1, fs_jpim1 ! vector opt.
	337	z3d(ji,jj,jk) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &
[12756]	338	& - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj)
[9403]	339	END DO
	340	END DO
	341	END DO
[10425]	342	CALL lbc_lnk( 'diawri', z3d, 'F', 1. )
[9403]	343	CALL iom_put( "relvor", z3d ) ! relative vorticity
	344
	345	DO jk = 1, jpkm1
	346	DO jj = 1, jpj
	347	DO ji = 1, jpi
	348	z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)
	349	END DO
	350	END DO
	351	END DO
	352	CALL iom_put( "absvor", z3d ) ! absolute vorticity
	353
	354	DO jk = 1, jpkm1
	355	DO jj = 1, jpjm1
	356	DO ji = 1, fs_jpim1 ! vector opt.
	357	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) &
	358	& + e3t_n(ji,jj ,jk)tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)tmask(ji+1,jj ,jk) )
	359	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	360	ELSE ; ze3 = 0._wp
	361	ENDIF
	362	z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)
	363	END DO
	364	END DO
	365	END DO
[10425]	366	CALL lbc_lnk( 'diawri', z3d, 'F', 1. )
[9403]	367	CALL iom_put( "potvor", z3d ) ! potential vorticity
	368
	369	ENDIF
	370
	371	!
	372	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
	373	z3d(:,:,jpk) = 0.e0
	374	z2d(:,:) = 0.e0
	375	DO jk = 1, jpkm1
	376	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
	377	z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
	378	END DO
	379	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
	380	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
	381	ENDIF
	382
	383	IF( iom_use("u_heattr") ) THEN
	384	z2d(:,:) = 0._wp
	385	DO jk = 1, jpkm1
	386	DO jj = 2, jpjm1
	387	DO ji = fs_2, fs_jpim1 ! vector opt.
	388	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
	389	END DO
	390	END DO
	391	END DO
[10425]	392	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9403]	393	CALL iom_put( "u_heattr", 0.5rcp z2d ) ! heat transport in i-direction
	394	ENDIF
	395
	396	IF( iom_use("u_salttr") ) THEN
	397	z2d(:,:) = 0.e0
	398	DO jk = 1, jpkm1
	399	DO jj = 2, jpjm1
	400	DO ji = fs_2, fs_jpim1 ! vector opt.
	401	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
	402	END DO
	403	END DO
	404	END DO
[10425]	405	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9403]	406	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
	407	ENDIF
	408
	409
	410	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
	411	z3d(:,:,jpk) = 0.e0
	412	DO jk = 1, jpkm1
	413	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
	414	END DO
	415	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
	416	ENDIF
	417
	418	IF( iom_use("v_heattr") ) THEN
	419	z2d(:,:) = 0.e0
	420	DO jk = 1, jpkm1
	421	DO jj = 2, jpjm1
	422	DO ji = fs_2, fs_jpim1 ! vector opt.
	423	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
	424	END DO
	425	END DO
	426	END DO
[10425]	427	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9403]	428	CALL iom_put( "v_heattr", 0.5rcp z2d ) ! heat transport in j-direction
	429	ENDIF
	430
	431	IF( iom_use("v_salttr") ) THEN
	432	z2d(:,:) = 0._wp
	433	DO jk = 1, jpkm1
	434	DO jj = 2, jpjm1
	435	DO ji = fs_2, fs_jpim1 ! vector opt.
	436	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
	437	END DO
	438	END DO
	439	END DO
[10425]	440	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9403]	441	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
	442	ENDIF
	443
	444	IF( iom_use("tosmint") ) THEN
	445	z2d(:,:) = 0._wp
	446	DO jk = 1, jpkm1
	447	DO jj = 2, jpjm1
	448	DO ji = fs_2, fs_jpim1 ! vector opt.
	449	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
	450	END DO
	451	END DO
	452	END DO
[10425]	453	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9403]	454	CALL iom_put( "tosmint", rau0 * z2d ) ! Vertical integral of temperature
	455	ENDIF
	456	IF( iom_use("somint") ) THEN
	457	z2d(:,:)=0._wp
	458	DO jk = 1, jpkm1
	459	DO jj = 2, jpjm1
	460	DO ji = fs_2, fs_jpim1 ! vector opt.
	461	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
	462	END DO
	463	END DO
	464	END DO
[10425]	465	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9403]	466	CALL iom_put( "somint", rau0 * z2d ) ! Vertical integral of salinity
	467	ENDIF
	468
	469	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
	470	!
	471
	472	IF (ln_dia25h) CALL dia_25h( kt ) ! 25h averaging
	473
	474	IF( ln_timing ) CALL timing_stop('dia_wri')
	475	!
	476	END SUBROUTINE dia_wri
	477
	478	#else
	479	!!----------------------------------------------------------------------
	480	!! Default option use IOIPSL library
	481	!!----------------------------------------------------------------------
[10425]	482
[9652]	483	INTEGER FUNCTION dia_wri_alloc()
	484	!!----------------------------------------------------------------------
	485	INTEGER, DIMENSION(2) :: ierr
	486	!!----------------------------------------------------------------------
	487	ierr = 0
	488	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	489	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	490	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
	491	!
	492	dia_wri_alloc = MAXVAL(ierr)
[10425]	493	CALL mpp_sum( 'diawri', dia_wri_alloc )
[9652]	494	!
	495	END FUNCTION dia_wri_alloc
[10425]	496
	497
[9403]	498	SUBROUTINE dia_wri( kt )
	499	!!---------------------------------------------------------------------
	500	!! * ROUTINE dia_wri *
	501	!!
	502	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	503	!! NETCDF format is used by default
	504	!!
	505	!! ** Method : At the beginning of the first time step (nit000),
	506	!! define all the NETCDF files and fields
	507	!! At each time step call histdef to compute the mean if ncessary
[11536]	508	!! Each nn_write time step, output the instantaneous or mean fields
[9403]	509	!!----------------------------------------------------------------------
	510	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	511	!
	512	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	513	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	514	INTEGER :: inum = 11 ! temporary logical unit
	515	INTEGER :: ji, jj, jk ! dummy loop indices
	516	INTEGER :: ierr ! error code return from allocation
	517	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
	518	INTEGER :: jn, ierror ! local integers
	519	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
	520	!
	521	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
	522	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
	523	!!----------------------------------------------------------------------
	524	!
	525	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[10425]	526	CALL dia_wri_state( 'output.init' )
[9403]	527	ninist = 0
	528	ENDIF
	529	!
[11536]	530	IF( nn_write == -1 ) RETURN ! we will never do any output
	531	!
	532	IF( ln_timing ) CALL timing_start('dia_wri')
	533	!
[9403]	534	! 0. Initialisation
	535	! -----------------
	536
	537	ll_print = .FALSE. ! local variable for debugging
	538	ll_print = ll_print .AND. lwp
	539
	540	! Define frequency of output and means
	541	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
	542	#if defined key_diainstant
[11536]	543	zsto = nn_write * rdt
[9403]	544	clop = "inst("//TRIM(clop)//")"
	545	#else
	546	zsto=rdt
	547	clop = "ave("//TRIM(clop)//")"
	548	#endif
[11536]	549	zout = nn_write * rdt
[9403]	550	zmax = ( nitend - nit000 + 1 ) * rdt
	551
	552	! Define indices of the horizontal output zoom and vertical limit storage
	553	iimi = 1 ; iima = jpi
	554	ijmi = 1 ; ijma = jpj
	555	ipk = jpk
	556
	557	! define time axis
	558	it = kt
	559	itmod = kt - nit000 + 1
	560
	561
	562	! 1. Define NETCDF files and fields at beginning of first time step
	563	! -----------------------------------------------------------------
	564
	565	IF( kt == nit000 ) THEN
	566
	567	! Define the NETCDF files (one per grid)
	568
	569	! Compute julian date from starting date of the run
	570	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	571	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
	572	IF(lwp)WRITE(numout,*)
	573	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	574	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	575	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	576	' limit storage in depth = ', ipk
	577
	578	! WRITE root name in date.file for use by postpro
	579	IF(lwp) THEN
[11536]	580	CALL dia_nam( clhstnam, nn_write,' ' )
[9403]	581	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
	582	WRITE(inum,*) clhstnam
	583	CLOSE(inum)
	584	ENDIF
	585
	586	! Define the T grid FILE ( nid_T )
	587
[11536]	588	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[9403]	589	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	590	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	591	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	592	& nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
	593	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
	594	& "m", ipk, gdept_1d, nz_T, "down" )
	595	! ! Index of ocean points
	596	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	597	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
	598	!
	599	IF( ln_icebergs ) THEN
	600	!
	601	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	602	!! that routine is called from nemogcm, so do it here immediately before its needed
	603	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	604	CALL mpp_sum( 'diawri', ierror )
[9403]	605	IF( ierror /= 0 ) THEN
	606	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	607	RETURN
	608	ENDIF
	609	!
	610	!! iceberg vertical coordinate is class number
	611	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	612	& "number", nclasses, class_num, nb_T )
	613	!
	614	!! each class just needs the surface index pattern
	615	ndim_bT = 3
	616	DO jn = 1,nclasses
	617	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	618	ENDDO
	619	!
	620	ENDIF
	621
	622	! Define the U grid FILE ( nid_U )
	623
[11536]	624	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[9403]	625	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	626	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	627	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	628	& nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
	629	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
	630	& "m", ipk, gdept_1d, nz_U, "down" )
	631	! ! Index of ocean points
	632	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	633	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	634
	635	! Define the V grid FILE ( nid_V )
	636
[11536]	637	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[9403]	638	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	639	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	640	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	641	& nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
	642	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
	643	& "m", ipk, gdept_1d, nz_V, "down" )
	644	! ! Index of ocean points
	645	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	646	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	647
	648	! Define the W grid FILE ( nid_W )
	649
[11536]	650	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[9403]	651	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	652	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	653	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	654	& nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
	655	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
	656	& "m", ipk, gdepw_1d, nz_W, "down" )
	657
	658
	659	! Declare all the output fields as NETCDF variables
	660
	661	! !!! nid_T : 3D
	662	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	663	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	664	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	665	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	666	IF( .NOT.ln_linssh ) THEN
	667	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	668	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	669	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	670	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	671	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	672	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	673	ENDIF
	674	! !!! nid_T : 2D
	675	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	676	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	677	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	678	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	679	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
	680	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	681	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
	682	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	683	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	684	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	685	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
	686	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	687	IF( ln_linssh ) THEN
	688	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
	689	& , "KgC/m2/s", & ! sosst_cd
	690	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	691	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
	692	& , "KgPSU/m2/s",& ! sosss_cd
	693	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	694	ENDIF
	695	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
	696	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	697	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	698	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	699	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	700	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	701	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	702	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	703	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
	704	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	705	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	706	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	707	!
	708	IF( ln_icebergs ) THEN
	709	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	710	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	711	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	712	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	713	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	714	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	715	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	716	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	717	IF( ln_bergdia ) THEN
	718	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	719	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	720	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	721	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	722	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	723	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	724	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	725	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	726	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	727	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	728	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	729	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	730	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	731	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	732	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	733	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	734	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	735	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	736	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	737	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	738	ENDIF
	739	ENDIF
	740
[11536]	741	IF( ln_ssr ) THEN
[9403]	742	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	743	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	744	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	745	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	746	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	747	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	748	ENDIF
[11536]	749
[9403]	750	clmx ="l_max(only(x))" ! max index on a period
	751	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	752	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
	753	#if defined key_diahth
	754	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	755	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	756	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	757	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	758	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	759	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	760	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
	761	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	762	#endif
	763
	764	CALL histend( nid_T, snc4chunks=snc4set )
	765
	766	! !!! nid_U : 3D
	767	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	768	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	769	IF( ln_wave .AND. ln_sdw) THEN
	770	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	771	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	772	ENDIF
	773	! !!! nid_U : 2D
	774	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
	775	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	776
	777	CALL histend( nid_U, snc4chunks=snc4set )
	778
	779	! !!! nid_V : 3D
	780	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	781	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	782	IF( ln_wave .AND. ln_sdw) THEN
	783	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	784	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	785	ENDIF
	786	! !!! nid_V : 2D
	787	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
	788	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	789
	790	CALL histend( nid_V, snc4chunks=snc4set )
	791
	792	! !!! nid_W : 3D
	793	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	794	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	795	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	796	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	797	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
	798	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	799
	800	IF( ln_zdfddm ) THEN
	801	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	802	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	803	ENDIF
	804
	805	IF( ln_wave .AND. ln_sdw) THEN
	806	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	807	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	808	ENDIF
	809	! !!! nid_W : 2D
	810	CALL histend( nid_W, snc4chunks=snc4set )
	811
	812	IF(lwp) WRITE(numout,*)
	813	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	814	IF(ll_print) CALL FLUSH(numout )
	815
	816	ENDIF
	817
	818	! 2. Start writing data
	819	! ---------------------
	820
	821	! ndex(1) est utilise ssi l'avant dernier argument est different de
	822	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	823	! donne le nombre d'elements, et ndex la liste des indices a sortir
	824
[11536]	825	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
[9403]	826	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	827	WRITE(numout,*) '~~~~~~ '
	828	ENDIF
	829
	830	IF( .NOT.ln_linssh ) THEN
	831	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	832	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	833	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	834	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
	835	ELSE
	836	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	837	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	838	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	839	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	840	ENDIF
	841	IF( .NOT.ln_linssh ) THEN
	842	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
	843	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	844	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
	845	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	846	ENDIF
	847	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
	848	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
	849	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
	850	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	851	! (includes virtual salt flux beneath ice
	852	! in linear free surface case)
	853	IF( ln_linssh ) THEN
	854	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
	855	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
	856	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
	857	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	858	ENDIF
	859	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
	860	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
	861	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
	862	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
	863	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
	864	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
	865	!
	866	IF( ln_icebergs ) THEN
	867	!
	868	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	869	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	870	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	871	!
	872	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	873	!
	874	IF( ln_bergdia ) THEN
	875	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	876	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	877	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	878	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	879	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	880	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	881	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	882	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	883	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	884	!
	885	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	886	ENDIF
	887	ENDIF
	888
[11536]	889	IF( ln_ssr ) THEN
[9403]	890	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	891	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[11536]	892	zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[9403]	893	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	894	ENDIF
	895	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	896	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
	897
	898	#if defined key_diahth
	899	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	900	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	901	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	902	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	903	#endif
	904
	905	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
	906	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
	907
	908	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
	909	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
	910
	911	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
	912	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
	913	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	914	IF( ln_zdfddm ) THEN
	915	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
	916	ENDIF
	917
	918	IF( ln_wave .AND. ln_sdw ) THEN
	919	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	920	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	921	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
	922	ENDIF
	923
	924	! 3. Close all files
	925	! ---------------------------------------
	926	IF( kt == nitend ) THEN
	927	CALL histclo( nid_T )
	928	CALL histclo( nid_U )
	929	CALL histclo( nid_V )
	930	CALL histclo( nid_W )
	931	ENDIF
	932	!
	933	IF( ln_timing ) CALL timing_stop('dia_wri')
	934	!
	935	END SUBROUTINE dia_wri
	936	#endif
	937
[10425]	938	SUBROUTINE dia_wri_state( cdfile_name )
[9403]	939	!!---------------------------------------------------------------------
	940	!! * ROUTINE dia_wri_state *
	941	!!
	942	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	943	!! the instantaneous ocean state and forcing fields.
	944	!! Used to find errors in the initial state or save the last
	945	!! ocean state in case of abnormal end of a simulation
	946	!!
	947	!! ** Method : NetCDF files using ioipsl
	948	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	949	!! File 'output.abort.nc' is created in case of abnormal job end
	950	!!----------------------------------------------------------------------
	951	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	952	!!
	953	INTEGER :: inum
[9403]	954	!!----------------------------------------------------------------------
	955	!
	956	IF(lwp) WRITE(numout,*)
	957	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	958	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
[10425]	959	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
[9403]	960
[9576]	961	#if defined key_si3
[10425]	962	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
[9403]	963	#else
[10425]	964	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
[9403]	965	#endif
	966
[10425]	967	CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) ) ! now temperature
	968	CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) ) ! now salinity
	969	CALL iom_rstput( 0, 0, inum, 'sossheig', sshn ) ! sea surface height
	970	CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity
	971	CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity
	972	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity
[9403]	973	IF( ALLOCATED(ahtu) ) THEN
[10425]	974	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	975	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[9403]	976	ENDIF
	977	IF( ALLOCATED(ahmt) ) THEN
[10425]	978	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	979	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[9403]	980	ENDIF
[10425]	981	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	982	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	983	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	984	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	985	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	986	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[9403]	987	IF( .NOT.ln_linssh ) THEN
[10425]	988	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
	989	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
	990	END IF
[9403]	991	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	992	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	993	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	994	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[9403]	995	ENDIF
[10425]	996
	997	#if defined key_si3
	998	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
	999	CALL ice_wri_state( inum )
[9403]	1000	ENDIF
	1001	#endif
[10425]	1002	!
	1003	CALL iom_close( inum )
[9403]	1004	!
	1005	END SUBROUTINE dia_wri_state
	1006
	1007	!!======================================================================
	1008	END MODULE diawri

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

Download in other formats: