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_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DIA/diawri.F90 @ 14644

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

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

Download in other formats: