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/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90 @ 11395

Last change on this file since 11395 was 11395, checked in by mathiot, 5 years ago
ENHANCE-02_ISF_nemo : Initial commit isf simplification (add ISF directory, moved isf routine in and split isf cavity and isf parametrisation, ...) (ticket #2142)
Property svn:keywords set to `Id`
File size: 50.5 KB

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

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

Download in other formats: