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/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DIA/diawri.F90 @ 13696

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

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

Download in other formats: