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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diawri.F90 in NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90 @ 12165

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

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

Download in other formats: