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 branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 11277

Last change on this file since 11277 was 11277, checked in by kingr, 5 years ago
Merged Juan's changes for running AMM15 woth wave coupling. Corrected minor logic error to allow AMM7-uncoupled to reproduce earlier results. Few line spacing changes to allow merging with OBS br and trunk rvn 5518.
File size: 59.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
	19	!!----------------------------------------------------------------------
[3]	20
	21	!!----------------------------------------------------------------------
[2528]	22	!! dia_wri : create the standart output files
	23	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
	24	!!----------------------------------------------------------------------
[3]	25	USE oce ! ocean dynamics and tracers
	26	USE dom_oce ! ocean space and time domain
[4292]	27	USE dynadv, ONLY: ln_dynadv_vec
[3]	28	USE zdf_oce ! ocean vertical physics
	29	USE ldftra_oce ! ocean active tracers: lateral physics
	30	USE ldfdyn_oce ! ocean dynamics: lateral physics
[3294]	31	USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
[3]	32	USE sol_oce ! solver variables
[888]	33	USE sbc_oce ! Surface boundary condition: ocean fields
	34	USE sbc_ice ! Surface boundary condition: ice fields
[3609]	35	USE icb_oce ! Icebergs
	36	USE icbdia ! Iceberg budgets
[888]	37	USE sbcssr ! restoring term toward SST/SSS climatology
[3]	38	USE phycst ! physical constants
	39	USE zdfmxl ! mixed layer
	40	USE dianam ! build name of file (routine)
	41	USE zdfddm ! vertical physics: double diffusion
	42	USE diahth ! thermocline diagnostics
[11277]	43	USE sbcwave ! wave parameters
[3]	44	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	45	USE in_out_manager ! I/O manager
[216]	46	USE diadimg ! dimg direct access file format output
[8059]	47	USE diatmb ! Top,middle,bottom output
	48	USE dia25h ! 25h Mean output
[8561]	49	USE diaopfoam ! Diaopfoam output
[1359]	50	USE iom
[460]	51	USE ioipsl
[5463]	52	USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities
[8059]	53	USE eosbn2 ! equation of state (eos_bn2 routine)
[5463]	54
[8059]	55
[1482]	56	#if defined key_lim2
	57	USE limwri_2
[4161]	58	#elif defined key_lim3
	59	USE limwri
[1482]	60	#endif
[2715]	61	USE lib_mpp ! MPP library
[3294]	62	USE timing ! preformance summary
	63	USE wrk_nemo ! working array
[2528]	64
[3]	65	IMPLICIT NONE
	66	PRIVATE
	67
[2528]	68	PUBLIC dia_wri ! routines called by step.F90
	69	PUBLIC dia_wri_state
[2715]	70	PUBLIC dia_wri_alloc ! Called by nemogcm module
[3]	71
[2528]	72	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[3609]	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
	77	INTEGER :: ndex(1) ! ???
[2715]	78	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
	79	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]	80	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	81
	82	!! * Substitutions
	83	# include "zdfddm_substitute.h90"
[1756]	84	# include "domzgr_substitute.h90"
	85	# include "vectopt_loop_substitute.h90"
[3]	86	!!----------------------------------------------------------------------
[2528]	87	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[5217]	88	!! $Id$
[2528]	89	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	90	!!----------------------------------------------------------------------
	91	CONTAINS
	92
[2715]	93	INTEGER FUNCTION dia_wri_alloc()
	94	!!----------------------------------------------------------------------
	95	INTEGER, DIMENSION(2) :: ierr
	96	!!----------------------------------------------------------------------
	97	ierr = 0
	98	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	99	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	100	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
	101	!
	102	dia_wri_alloc = MAXVAL(ierr)
	103	IF( lk_mpp ) CALL mpp_sum( dia_wri_alloc )
	104	!
	105	END FUNCTION dia_wri_alloc
	106
[216]	107	#if defined key_dimgout
[3]	108	!!----------------------------------------------------------------------
[2528]	109	!! 'key_dimgout' DIMG output file
[3]	110	!!----------------------------------------------------------------------
	111	# include "diawri_dimg.h90"
	112
	113	#else
	114	!!----------------------------------------------------------------------
	115	!! Default option NetCDF output file
	116	!!----------------------------------------------------------------------
[2528]	117	# if defined key_iomput
[3]	118	!!----------------------------------------------------------------------
[2528]	119	!! 'key_iomput' use IOM library
	120	!!----------------------------------------------------------------------
[2715]	121
[1561]	122	SUBROUTINE dia_wri( kt )
[1482]	123	!!---------------------------------------------------------------------
	124	!! * ROUTINE dia_wri *
	125	!!
	126	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	127	!! NETCDF format is used by default
	128	!!
	129	!! ** Method : use iom_put
	130	!!----------------------------------------------------------------------
[1756]	131	!!
[1482]	132	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1756]	133	!!
	134	INTEGER :: ji, jj, jk ! dummy loop indices
[5460]	135	INTEGER :: jkbot !
[1756]	136	REAL(wp) :: zztmp, zztmpx, zztmpy !
[3294]	137	!!
[4761]	138	REAL(wp), POINTER, DIMENSION(:,:) :: z2d ! 2D workspace
[3294]	139	REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d ! 3D workspace
[8059]	140	REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhd , zrhop ! 3D workspace
[1482]	141	!!----------------------------------------------------------------------
	142	!
[3294]	143	IF( nn_timing == 1 ) CALL timing_start('dia_wri')
	144	!
[4990]	145	CALL wrk_alloc( jpi , jpj , z2d )
[3294]	146	CALL wrk_alloc( jpi , jpj, jpk , z3d )
[8059]	147	CALL wrk_alloc( jpi , jpj, jpk , zrhd , zrhop )
[2715]	148	!
[1482]	149	! Output the initial state and forcings
	150	IF( ninist == 1 ) THEN
	151	CALL dia_wri_state( 'output.init', kt )
	152	ninist = 0
	153	ENDIF
[3]	154
[5107]	155	IF( .NOT.lk_vvl ) THEN
	156	CALL iom_put( "e3t" , fse3t_n(:,:,:) )
	157	CALL iom_put( "e3u" , fse3u_n(:,:,:) )
	158	CALL iom_put( "e3v" , fse3v_n(:,:,:) )
	159	CALL iom_put( "e3w" , fse3w_n(:,:,:) )
	160	ENDIF
[5461]	161
	162	CALL iom_put( "ssh" , sshn ) ! sea surface height
	163	if( iom_use('ssh2') ) CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height
[5107]	164
	165	CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature
	166	CALL iom_put( "sst", tsn(:,:,1,jp_tem) ) ! surface temperature
	167	IF ( iom_use("sbt") ) THEN
[4990]	168	DO jj = 1, jpj
	169	DO ji = 1, jpi
[5460]	170	jkbot = mbkt(ji,jj)
	171	z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
[4990]	172	END DO
[5107]	173	END DO
	174	CALL iom_put( "sbt", z2d ) ! bottom temperature
	175	ENDIF
	176
	177	CALL iom_put( "soce", tsn(:,:,:,jp_sal) ) ! 3D salinity
	178	CALL iom_put( "sss", tsn(:,:,1,jp_sal) ) ! surface salinity
	179	IF ( iom_use("sbs") ) THEN
[4990]	180	DO jj = 1, jpj
	181	DO ji = 1, jpi
[5460]	182	jkbot = mbkt(ji,jj)
	183	z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
[4990]	184	END DO
[5107]	185	END DO
	186	CALL iom_put( "sbs", z2d ) ! bottom salinity
	187	ENDIF
[5463]	188
	189	IF ( iom_use("taubot") ) THEN ! bottom stress
	190	z2d(:,:) = 0._wp
	191	DO jj = 2, jpjm1
	192	DO ji = fs_2, fs_jpim1 ! vector opt.
	193	zztmpx = ( bfrua(ji ,jj) * un(ji ,jj,mbku(ji ,jj)) &
	194	& + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj)) )
	195	zztmpy = ( bfrva(ji, jj) * vn(ji,jj ,mbkv(ji,jj )) &
	196	& + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1)) )
	197	z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)
	198	!
	199	ENDDO
	200	ENDDO
	201	CALL lbc_lnk( z2d, 'T', 1. )
	202	CALL iom_put( "taubot", z2d )
	203	ENDIF
[5107]	204
	205	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
	206	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
	207	IF ( iom_use("sbu") ) THEN
[4990]	208	DO jj = 1, jpj
	209	DO ji = 1, jpi
[5460]	210	jkbot = mbku(ji,jj)
	211	z2d(ji,jj) = un(ji,jj,jkbot)
[4990]	212	END DO
[5107]	213	END DO
	214	CALL iom_put( "sbu", z2d ) ! bottom i-current
	215	ENDIF
[5463]	216	#if defined key_dynspg_ts
	217	CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current
	218	#else
	219	CALL iom_put( "ubar", un_b(:,:) ) ! barotropic i-current
	220	#endif
[5107]	221
	222	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
	223	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
	224	IF ( iom_use("sbv") ) THEN
[4990]	225	DO jj = 1, jpj
	226	DO ji = 1, jpi
[5460]	227	jkbot = mbkv(ji,jj)
	228	z2d(ji,jj) = vn(ji,jj,jkbot)
[4990]	229	END DO
[5107]	230	END DO
	231	CALL iom_put( "sbv", z2d ) ! bottom j-current
[4990]	232	ENDIF
[5463]	233	#if defined key_dynspg_ts
	234	CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic j-current
	235	#else
	236	CALL iom_put( "vbar", vn_b(:,:) ) ! barotropic j-current
	237	#endif
[1482]	238
[5461]	239	CALL iom_put( "woce", wn ) ! vertical velocity
	240	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
	241	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
	242	z2d(:,:) = rau0 * e12t(:,:)
	243	DO jk = 1, jpk
	244	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
	245	END DO
	246	CALL iom_put( "w_masstr" , z3d )
	247	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
	248	ENDIF
	249
[5107]	250	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	251	CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef.
	252	CALL iom_put( "avs" , fsavs(:,:,:) ) ! S vert. eddy diff. coef. (useful only with key_zdfddm)
	253
	254	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
[4990]	255	DO jj = 2, jpjm1 ! sst gradient
	256	DO ji = fs_2, fs_jpim1 ! vector opt.
	257	zztmp = tsn(ji,jj,1,jp_tem)
	258	zztmpx = ( tsn(ji+1,jj ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,1,jp_tem) ) / e1u(ji-1,jj )
	259	zztmpy = ( tsn(ji ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,1,jp_tem) ) / e2v(ji ,jj-1)
	260	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	261	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
	262	END DO
[1756]	263	END DO
[4990]	264	CALL lbc_lnk( z2d, 'T', 1. )
	265	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
	266	z2d(:,:) = SQRT( z2d(:,:) )
	267	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
	268	ENDIF
	269
[4761]	270	! clem: heat and salt content
[4990]	271	IF( iom_use("heatc") ) THEN
	272	z2d(:,:) = 0._wp
	273	DO jk = 1, jpkm1
[5107]	274	DO jj = 1, jpj
	275	DO ji = 1, jpi
[4990]	276	z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
	277	END DO
[4761]	278	END DO
	279	END DO
[4990]	280	CALL iom_put( "heatc", (rau0 * rcp) * z2d ) ! vertically integrated heat content (J/m2)
	281	ENDIF
	282
	283	IF( iom_use("saltc") ) THEN
	284	z2d(:,:) = 0._wp
	285	DO jk = 1, jpkm1
[5107]	286	DO jj = 1, jpj
	287	DO ji = 1, jpi
[4990]	288	z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
	289	END DO
	290	END DO
	291	END DO
	292	CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
	293	ENDIF
[4840]	294	!
[4990]	295	IF ( iom_use("eken") ) THEN
	296	rke(:,:,jk) = 0._wp ! kinetic energy
	297	DO jk = 1, jpkm1
	298	DO jj = 2, jpjm1
	299	DO ji = fs_2, fs_jpim1 ! vector opt.
	300	zztmp = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
	301	zztmpx = 0.5 * ( un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) &
	302	& + un(ji ,jj,jk) * un(ji ,jj,jk) * e2u(ji ,jj) * fse3u(ji ,jj,jk) ) &
	303	& * zztmp
	304	!
	305	zztmpy = 0.5 * ( vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) &
	306	& + vn(ji,jj ,jk) * vn(ji,jj ,jk) * e1v(ji,jj ) * fse3v(ji,jj ,jk) ) &
	307	& * zztmp
	308	!
	309	rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
	310	!
	311	ENDDO
[4840]	312	ENDDO
	313	ENDDO
[4990]	314	CALL lbc_lnk( rke, 'T', 1. )
	315	CALL iom_put( "eken", rke )
	316	ENDIF
	317
	318	IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
[1756]	319	z3d(:,:,jpk) = 0.e0
	320	DO jk = 1, jpkm1
[4761]	321	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
[1756]	322	END DO
	323	CALL iom_put( "u_masstr", z3d ) ! mass transport in i-direction
[4990]	324	ENDIF
	325
	326	IF( iom_use("u_heattr") ) THEN
	327	z2d(:,:) = 0.e0
	328	DO jk = 1, jpkm1
	329	DO jj = 2, jpjm1
	330	DO ji = fs_2, fs_jpim1 ! vector opt.
	331	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
	332	END DO
	333	END DO
	334	END DO
	335	CALL lbc_lnk( z2d, 'U', -1. )
	336	CALL iom_put( "u_heattr", (0.5 * rcp) * z2d ) ! heat transport in i-direction
	337	ENDIF
[4761]	338
[4990]	339	IF( iom_use("u_salttr") ) THEN
[1756]	340	z2d(:,:) = 0.e0
	341	DO jk = 1, jpkm1
	342	DO jj = 2, jpjm1
	343	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	344	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
[1756]	345	END DO
	346	END DO
	347	END DO
	348	CALL lbc_lnk( z2d, 'U', -1. )
[4990]	349	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
	350	ENDIF
[4761]	351
[4990]	352
	353	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
[4761]	354	z3d(:,:,jpk) = 0.e0
[1756]	355	DO jk = 1, jpkm1
[4761]	356	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
[1756]	357	END DO
	358	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
[4990]	359	ENDIF
	360
	361	IF( iom_use("v_heattr") ) THEN
	362	z2d(:,:) = 0.e0
	363	DO jk = 1, jpkm1
	364	DO jj = 2, jpjm1
	365	DO ji = fs_2, fs_jpim1 ! vector opt.
	366	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
	367	END DO
	368	END DO
	369	END DO
	370	CALL lbc_lnk( z2d, 'V', -1. )
	371	CALL iom_put( "v_heattr", (0.5 * rcp) * z2d ) ! heat transport in j-direction
	372	ENDIF
[4761]	373
[4990]	374	IF( iom_use("v_salttr") ) THEN
[1756]	375	z2d(:,:) = 0.e0
	376	DO jk = 1, jpkm1
	377	DO jj = 2, jpjm1
	378	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	379	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
[1756]	380	END DO
	381	END DO
	382	END DO
	383	CALL lbc_lnk( z2d, 'V', -1. )
[4990]	384	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
[1756]	385	ENDIF
[8059]	386
	387	IF( iom_use("rhop") ) THEN
	388	CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) ) ! now in situ and potential density
	389	zrhop(:,:,jpk) = 0._wp
	390	CALL iom_put( 'rhop', zrhop )
	391	ENDIF
	392
[2528]	393	!
[4990]	394	CALL wrk_dealloc( jpi , jpj , z2d )
[3294]	395	CALL wrk_dealloc( jpi , jpj, jpk , z3d )
[8059]	396	CALL wrk_dealloc( jpi , jpj, jpk , zrhd , zrhop )
[2715]	397	!
[8059]	398	! If we want tmb values
	399
	400	IF (ln_diatmb) THEN
	401	CALL dia_tmb
	402	ENDIF
	403	IF (ln_dia25h) THEN
	404	CALL dia_25h( kt )
	405	ENDIF
[8561]	406	IF (ln_diaopfoam) THEN
	407	CALL dia_diaopfoam
	408	ENDIF
[8059]	409	!
[3294]	410	IF( nn_timing == 1 ) CALL timing_stop('dia_wri')
	411	!
[1482]	412	END SUBROUTINE dia_wri
	413
	414	#else
[2528]	415	!!----------------------------------------------------------------------
	416	!! Default option use IOIPSL library
	417	!!----------------------------------------------------------------------
	418
[1561]	419	SUBROUTINE dia_wri( kt )
[3]	420	!!---------------------------------------------------------------------
	421	!! * ROUTINE dia_wri *
	422	!!
	423	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	424	!! NETCDF format is used by default
	425	!!
	426	!! ** Method : At the beginning of the first time step (nit000),
	427	!! define all the NETCDF files and fields
	428	!! At each time step call histdef to compute the mean if ncessary
	429	!! Each nwrite time step, output the instantaneous or mean fields
	430	!!----------------------------------------------------------------------
[2715]	431	!!
[3]	432	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2528]	433	!!
	434	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	435	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	436	INTEGER :: inum = 11 ! temporary logical unit
[3294]	437	INTEGER :: ji, jj, jk ! dummy loop indices
	438	INTEGER :: ierr ! error code return from allocation
[2528]	439	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[3609]	440	INTEGER :: jn, ierror ! local integers
[2528]	441	REAL(wp) :: zsto, zout, zmax, zjulian, zdt ! local scalars
[3294]	442	!!
	443	REAL(wp), POINTER, DIMENSION(:,:) :: zw2d ! 2D workspace
	444	REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d ! 3D workspace
[3]	445	!!----------------------------------------------------------------------
[3294]	446	!
	447	IF( nn_timing == 1 ) CALL timing_start('dia_wri')
[1482]	448	!
[3294]	449	CALL wrk_alloc( jpi , jpj , zw2d )
[4292]	450	IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_alloc( jpi , jpj , jpk , zw3d )
[2715]	451	!
[1482]	452	! Output the initial state and forcings
	453	IF( ninist == 1 ) THEN
	454	CALL dia_wri_state( 'output.init', kt )
	455	ninist = 0
	456	ENDIF
	457	!
[3]	458	! 0. Initialisation
	459	! -----------------
[632]	460
[3]	461	! local variable for debugging
	462	ll_print = .FALSE.
	463	ll_print = ll_print .AND. lwp
	464
	465	! Define frequency of output and means
	466	zdt = rdt
	467	IF( nacc == 1 ) zdt = rdtmin
[5565]	468	clop = "x" ! no use of the mask value (require less cpu time, and otherwise the model crashes)
[3]	469	#if defined key_diainstant
	470	zsto = nwrite * zdt
[1312]	471	clop = "inst("//TRIM(clop)//")"
[3]	472	#else
	473	zsto=zdt
[1312]	474	clop = "ave("//TRIM(clop)//")"
[3]	475	#endif
	476	zout = nwrite * zdt
	477	zmax = ( nitend - nit000 + 1 ) * zdt
	478
	479	! Define indices of the horizontal output zoom and vertical limit storage
	480	iimi = 1 ; iima = jpi
	481	ijmi = 1 ; ijma = jpj
	482	ipk = jpk
	483
	484	! define time axis
[1334]	485	it = kt
	486	itmod = kt - nit000 + 1
[3]	487
	488
	489	! 1. Define NETCDF files and fields at beginning of first time step
	490	! -----------------------------------------------------------------
	491
	492	IF( kt == nit000 ) THEN
	493
	494	! Define the NETCDF files (one per grid)
[632]	495
[3]	496	! Compute julian date from starting date of the run
[1309]	497	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	498	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	499	IF(lwp)WRITE(numout,*)
	500	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	501	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	502	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	503	' limit storage in depth = ', ipk
	504
	505	! WRITE root name in date.file for use by postpro
[1581]	506	IF(lwp) THEN
[895]	507	CALL dia_nam( clhstnam, nwrite,' ' )
[1581]	508	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	509	WRITE(inum,*) clhstnam
	510	CLOSE(inum)
	511	ENDIF
[632]	512
[3]	513	! Define the T grid FILE ( nid_T )
[632]	514
[3]	515	CALL dia_nam( clhstnam, nwrite, 'grid_T' )
	516	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	517	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	518	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	519	& nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	520	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	521	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	522	! ! Index of ocean points
	523	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	524	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	525	!
	526	IF( ln_icebergs ) THEN
	527	!
	528	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	529	!! that routine is called from nemogcm, so do it here immediately before its needed
	530	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
	531	IF( lk_mpp ) CALL mpp_sum( ierror )
	532	IF( ierror /= 0 ) THEN
	533	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	534	RETURN
	535	ENDIF
	536	!
	537	!! iceberg vertical coordinate is class number
	538	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	539	& "number", nclasses, class_num, nb_T )
	540	!
	541	!! each class just needs the surface index pattern
	542	ndim_bT = 3
	543	DO jn = 1,nclasses
	544	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	545	ENDDO
	546	!
	547	ENDIF
[3]	548
	549	! Define the U grid FILE ( nid_U )
	550
	551	CALL dia_nam( clhstnam, nwrite, 'grid_U' )
	552	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	553	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	554	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	555	& nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	556	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	557	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	558	! ! Index of ocean points
	559	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	560	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	561
	562	! Define the V grid FILE ( nid_V )
	563
	564	CALL dia_nam( clhstnam, nwrite, 'grid_V' ) ! filename
	565	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	566	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	567	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	568	& nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	569	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	570	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	571	! ! Index of ocean points
	572	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	573	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	574
	575	! Define the W grid FILE ( nid_W )
	576
	577	CALL dia_nam( clhstnam, nwrite, 'grid_W' ) ! filename
	578	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	579	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	580	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	581	& nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	582	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	583	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	584
[632]	585
[3]	586	! Declare all the output fields as NETCDF variables
	587
	588	! !!! nid_T : 3D
	589	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	590	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	591	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	592	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[4292]	593	IF( lk_vvl ) THEN
	594	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	595	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	596	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	597	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	598	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	599	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	600	ENDIF
[3]	601	! !!! nid_T : 2D
	602	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	603	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	604	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	605	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	606	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	607	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	608	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	609	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	610	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	611	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	612	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	613	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4292]	614	IF( .NOT. lk_vvl ) THEN
	615	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
[3625]	616	& , "KgC/m2/s", & ! sosst_cd
[4292]	617	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	618	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
[3625]	619	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	620	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	621	ENDIF
[888]	622	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	623	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	624	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	625	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	626	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	627	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	628	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	629	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	630	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	631	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	632	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	633	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3609]	634	!
	635	IF( ln_icebergs ) THEN
	636	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	637	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	638	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	639	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	640	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	641	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	642	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	643	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	644	IF( ln_bergdia ) THEN
	645	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	646	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	647	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	648	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	649	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	650	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	651	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	652	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	653	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	654	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	655	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	656	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	657	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	658	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	659	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	660	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	661	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	662	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	663	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	664	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	665	ENDIF
	666	ENDIF
	667
[5407]	668	IF( .NOT. ln_cpl ) THEN
[4990]	669	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	670	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	671	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	672	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	673	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	674	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	675	ENDIF
[3]	676
[5407]	677	IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]	678	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	679	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	680	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	681	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	682	CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping" , "Kg/m2/s", & ! erp * sn
	683	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	684	ENDIF
	685
[3]	686	clmx ="l_max(only(x))" ! max index on a period
	687	CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	688	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
	689	#if defined key_diahth
	690	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	691	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	692	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	693	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	694	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	695	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	696	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "W" , & ! htc3
	697	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	698	#endif
	699
[5407]	700	IF( ln_cpl .AND. nn_ice == 2 ) THEN
[4990]	701	CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature" , "K" , & ! tn_ice
	702	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	703	CALL histdef( nid_T,"soicealb" , "Ice Albedo" , "[0,1]" , & ! alb_ice
	704	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	705	ENDIF
[3]	706
[2528]	707	CALL histend( nid_T, snc4chunks=snc4set )
[3]	708
	709	! !!! nid_U : 3D
	710	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	711	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[11277]	712	IF( ln_wave .AND. ln_sdw) THEN
	713	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	714	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	715	ENDIF
[3294]	716	IF( ln_traldf_gdia ) THEN
	717	CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current" , "m/s" , & ! u_eiv
	718	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	719	ELSE
[3]	720	#if defined key_diaeiv
[3294]	721	CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current" , "m/s" , & ! u_eiv
[3]	722	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	723	#endif
[3294]	724	END IF
[3]	725	! !!! nid_U : 2D
[888]	726	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	727	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	728
[2528]	729	CALL histend( nid_U, snc4chunks=snc4set )
[3]	730
	731	! !!! nid_V : 3D
	732	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	733	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[11277]	734	IF( ln_wave .AND. ln_sdw) THEN
	735	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	736	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	737	ENDIF
[3294]	738	IF( ln_traldf_gdia ) THEN
	739	CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current" , "m/s" , & ! v_eiv
	740	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	741	ELSE
[3]	742	#if defined key_diaeiv
[3294]	743	CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current" , "m/s" , & ! v_eiv
[3]	744	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	745	#endif
[3294]	746	END IF
[3]	747	! !!! nid_V : 2D
[888]	748	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	749	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	750
[2528]	751	CALL histend( nid_V, snc4chunks=snc4set )
[3]	752
	753	! !!! nid_W : 3D
	754	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	755	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[3294]	756	IF( ln_traldf_gdia ) THEN
	757	CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity" , "m/s" , & ! w_eiv
	758	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	759	ELSE
[3]	760	#if defined key_diaeiv
[3294]	761	CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity" , "m/s" , & ! w_eiv
	762	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[3]	763	#endif
[3294]	764	END IF
[3]	765	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	766	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[255]	767	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avmu
	768	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	769
[3]	770	IF( lk_zdfddm ) THEN
	771	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	772	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	773	ENDIF
[11277]	774
	775	IF( ln_wave .AND. ln_sdw) THEN
	776	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	777	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	778	ENDIF
[3]	779	! !!! nid_W : 2D
	780	#if defined key_traldf_c2d
	781	CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity" , "m2/s" , & ! ahtw
	782	& jpi, jpj, nh_W, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
[23]	783	# if defined key_traldf_eiv
[3]	784	CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s", & ! aeiw
	785	& jpi, jpj, nh_W, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
[23]	786	# endif
[3]	787	#endif
	788
[2528]	789	CALL histend( nid_W, snc4chunks=snc4set )
[3]	790
	791	IF(lwp) WRITE(numout,*)
	792	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	793	IF(ll_print) CALL FLUSH(numout )
	794
	795	ENDIF
	796
	797	! 2. Start writing data
	798	! ---------------------
	799
[4292]	800	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	801	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	802	! donne le nombre d'elements, et ndex la liste des indices a sortir
	803
[1334]	804	IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
[3]	805	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	806	WRITE(numout,*) '~~~~~~ '
	807	ENDIF
	808
[4292]	809	IF( lk_vvl ) THEN
[4492]	810	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	811	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	812	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	813	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	814	ELSE
	815	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	816	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	817	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	818	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	819	ENDIF
	820	IF( lk_vvl ) THEN
	821	zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
	822	CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	823	CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
	824	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	825	ENDIF
[3]	826	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
[2528]	827	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	828	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	829	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	830	! (includes virtual salt flux beneath ice
	831	! in linear free surface case)
[4292]	832	IF( .NOT. lk_vvl ) THEN
	833	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
	834	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
	835	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
	836	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	837	ENDIF
[888]	838	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	839	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	840	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	841	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	842	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	843	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[3609]	844	!
	845	IF( ln_icebergs ) THEN
	846	!
	847	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	848	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	849	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	850	!
	851	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	852	!
	853	IF( ln_bergdia ) THEN
	854	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	855	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	856	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	857	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	858	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	859	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	860	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	861	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	862	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	863	!
	864	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	865	ENDIF
	866	ENDIF
	867
[5407]	868	IF( .NOT. ln_cpl ) THEN
[4990]	869	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	870	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[3294]	871	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	872	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	873	ENDIF
[5407]	874	IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]	875	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	876	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
	877	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
	878	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	879	ENDIF
	880	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	881	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	882
	883	#if defined key_diahth
	884	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	885	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	886	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	887	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	888	#endif
[888]	889
[5407]	890	IF( ln_cpl .AND. nn_ice == 2 ) THEN
[4990]	891	CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT ) ! surf. ice temperature
	892	CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT ) ! ice albedo
	893	ENDIF
	894
[3]	895	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
[3294]	896	IF( ln_traldf_gdia ) THEN
	897	IF (.not. ALLOCATED(psix_eiv))THEN
	898	ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
	899	IF( lk_mpp ) CALL mpp_sum ( ierr )
	900	IF( ierr > 0 ) CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
	901	psix_eiv(:,:,:) = 0.0_wp
	902	psiy_eiv(:,:,:) = 0.0_wp
	903	ENDIF
	904	DO jk=1,jpkm1
	905	zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk) ! u_eiv = -dpsix/dz
	906	END DO
	907	zw3d(:,:,jpk) = 0._wp
	908	CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U ) ! i-eiv current
	909	ELSE
[3]	910	#if defined key_diaeiv
[3294]	911	CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U ) ! i-eiv current
[3]	912	#endif
[3294]	913	ENDIF
[888]	914	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	915
	916	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
[3294]	917	IF( ln_traldf_gdia ) THEN
	918	DO jk=1,jpk-1
	919	zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk) ! v_eiv = -dpsiy/dz
	920	END DO
	921	zw3d(:,:,jpk) = 0._wp
	922	CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V ) ! j-eiv current
	923	ELSE
[3]	924	#if defined key_diaeiv
[3294]	925	CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V ) ! j-eiv current
[3]	926	#endif
[3294]	927	ENDIF
[888]	928	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	929
	930	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
[3294]	931	IF( ln_traldf_gdia ) THEN
	932	DO jk=1,jpk-1
	933	DO jj = 2, jpjm1
	934	DO ji = fs_2, fs_jpim1 ! vector opt.
	935	zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
	936	& (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
	937	END DO
	938	END DO
	939	END DO
	940	zw3d(:,:,jpk) = 0._wp
	941	CALL histwrite( nid_W, "voveeivw", it, zw3d , ndim_T, ndex_T ) ! vert. eiv current
	942	ELSE
[3]	943	# if defined key_diaeiv
[3294]	944	CALL histwrite( nid_W, "voveeivw", it, w_eiv , ndim_T, ndex_T ) ! vert. eiv current
[3]	945	# endif
[3294]	946	ENDIF
[3]	947	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[255]	948	CALL histwrite( nid_W, "votkeavm", it, avmu , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
[3]	949	IF( lk_zdfddm ) THEN
	950	CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T ) ! S vert. eddy diff. coef.
	951	ENDIF
	952	#if defined key_traldf_c2d
	953	CALL histwrite( nid_W, "soleahtw", it, ahtw , ndim_hT, ndex_hT ) ! lateral eddy diff. coef.
[23]	954	# if defined key_traldf_eiv
[3]	955	CALL histwrite( nid_W, "soleaeiw", it, aeiw , ndim_hT, ndex_hT ) ! EIV coefficient at w-point
[23]	956	# endif
[3]	957	#endif
	958
[11277]	959	IF( ln_wave .AND. ln_sdw ) THEN
	960	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	961	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	962	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
	963	ENDIF
	964
[1318]	965	! 3. Close all files
	966	! ---------------------------------------
[1561]	967	IF( kt == nitend ) THEN
[3]	968	CALL histclo( nid_T )
	969	CALL histclo( nid_U )
	970	CALL histclo( nid_V )
	971	CALL histclo( nid_W )
	972	ENDIF
[2528]	973	!
[3294]	974	CALL wrk_dealloc( jpi , jpj , zw2d )
[4292]	975	IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_dealloc( jpi , jpj , jpk , zw3d )
[2715]	976	!
[3294]	977	IF( nn_timing == 1 ) CALL timing_stop('dia_wri')
	978	!
[3]	979	END SUBROUTINE dia_wri
[1482]	980	# endif
[3]	981
[1567]	982	#endif
	983
[1334]	984	SUBROUTINE dia_wri_state( cdfile_name, kt )
[3]	985	!!---------------------------------------------------------------------
	986	!! * ROUTINE dia_wri_state *
	987	!!
	988	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	989	!! the instantaneous ocean state and forcing fields.
	990	!! Used to find errors in the initial state or save the last
	991	!! ocean state in case of abnormal end of a simulation
	992	!!
	993	!! ** Method : NetCDF files using ioipsl
	994	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	995	!! File 'output.abort.nc' is created in case of abnormal job end
	996	!!----------------------------------------------------------------------
[1334]	997	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
	998	INTEGER , INTENT( in ) :: kt ! ocean time-step index
[2528]	999	!!
[648]	1000	CHARACTER (len=32) :: clname
[3]	1001	CHARACTER (len=40) :: clop
[2528]	1002	INTEGER :: id_i , nz_i, nh_i
	1003	INTEGER, DIMENSION(1) :: idex ! local workspace
	1004	REAL(wp) :: zsto, zout, zmax, zjulian, zdt
[3]	1005	!!----------------------------------------------------------------------
[3294]	1006	!
[3625]	1007	! IF( nn_timing == 1 ) CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
[3]	1008
	1009	! 0. Initialisation
	1010	! -----------------
[632]	1011
[648]	1012	! Define name, frequency of output and means
	1013	clname = cdfile_name
[1792]	1014	IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
[3]	1015	zdt = rdt
	1016	zsto = rdt
	1017	clop = "inst(x)" ! no use of the mask value (require less cpu time)
	1018	zout = rdt
	1019	zmax = ( nitend - nit000 + 1 ) * zdt
	1020
[648]	1021	IF(lwp) WRITE(numout,*)
	1022	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	1023	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
	1024	IF(lwp) WRITE(numout,*) ' and named :', clname, '.nc'
	1025
	1026
[3]	1027	! 1. Define NETCDF files and fields at beginning of first time step
	1028	! -----------------------------------------------------------------
	1029
	1030	! Compute julian date from starting date of the run
[1310]	1031	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) ! time axis
	1032	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[648]	1033	CALL histbeg( clname, jpi, glamt, jpj, gphit, &
[2528]	1034	1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
[3]	1035	CALL histvert( id_i, "deptht", "Vertical T levels", & ! Vertical grid : gdept
[4292]	1036	"m", jpk, gdept_1d, nz_i, "down")
[3]	1037
	1038	! Declare all the output fields as NetCDF variables
	1039
	1040	CALL histdef( id_i, "vosaline", "Salinity" , "PSU" , & ! salinity
	1041	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1042	CALL histdef( id_i, "votemper", "Temperature" , "C" , & ! temperature
	1043	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[359]	1044	CALL histdef( id_i, "sossheig", "Sea Surface Height" , "m" , & ! ssh
	1045	& jpi, jpj, nh_i, 1 , 1, 1 , nz_i, 32, clop, zsto, zout )
[3]	1046	CALL histdef( id_i, "vozocrtx", "Zonal Current" , "m/s" , & ! zonal current
	1047	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1048	CALL histdef( id_i, "vomecrty", "Meridional Current" , "m/s" , & ! meridonal current
	1049	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1050	CALL histdef( id_i, "vovecrtz", "Vertical Velocity" , "m/s" , & ! vertical current
	1051	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1052	CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S", & ! net freshwater
	1053	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1054	CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2" , & ! net heat flux
	1055	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1056	CALL histdef( id_i, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! solar flux
	1057	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	1058	CALL histdef( id_i, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	1059	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1060	CALL histdef( id_i, "sozotaux", "Zonal Wind Stress" , "N/m2" , & ! i-wind stress
	1061	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1062	CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2" , & ! j-wind stress
	1063	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4292]	1064	IF( lk_vvl ) THEN
	1065	CALL histdef( id_i, "vovvldep", "T point depth" , "m" , & ! t-point depth
	1066	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[6204]	1067	CALL histdef( id_i, "vovvle3t", "T point thickness" , "m" , & ! t-point depth
	1068	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[4292]	1069	END IF
[11277]	1070	!
	1071	IF( ln_wave .AND. ln_sdw ) THEN
	1072	CALL histdef( id_i, "sdzocrtx", "Stokes Drift Zonal", "m/s" , & ! StokesDrift zonal current
	1073	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1074	CALL histdef( id_i, "sdmecrty", "Stokes Drift Merid", "m/s" , & ! StokesDrift meridonal current
	1075	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1076	CALL histdef( id_i, "sdvecrtz", "Stokes Drift Vert", "m/s" , & ! StokesDrift vertical current
	1077	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1078	ENDIF
[3]	1079
[1482]	1080	#if defined key_lim2
	1081	CALL lim_wri_state_2( kt, id_i, nh_i )
[4161]	1082	#elif defined key_lim3
	1083	CALL lim_wri_state( kt, id_i, nh_i )
[1482]	1084	#else
[2528]	1085	CALL histend( id_i, snc4chunks=snc4set )
[1482]	1086	#endif
[3]	1087
	1088	! 2. Start writing data
	1089	! ---------------------
	1090	! idex(1) est utilise ssi l'avant dernier argument est diffferent de
	1091	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	1092	! donne le nombre d'elements, et idex la liste des indices a sortir
	1093	idex(1) = 1 ! init to avoid compil warning
[632]	1094
[3]	1095	! Write all fields on T grid
[11277]	1096	IF( ln_wave .AND. ln_sdw ) THEN
	1097	CALL histwrite( id_i, "sdzocrtx", kt, usd , jpijpjjpk, idex) ! now StokesDrift i-velocity
	1098	CALL histwrite( id_i, "sdmecrty", kt, vsd , jpijpjjpk, idex) ! now StokesDrift j-velocity
	1099	CALL histwrite( id_i, "sdvecrtz", kt, wsd , jpijpjjpk, idex) ! now StokesDrift k-velocity
	1100	ENDIF
[3294]	1101	CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpijpjjpk, idex ) ! now temperature
	1102	CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpijpjjpk, idex ) ! now salinity
	1103	CALL histwrite( id_i, "sossheig", kt, sshn , jpi*jpj , idex ) ! sea surface height
	1104	CALL histwrite( id_i, "vozocrtx", kt, un , jpijpjjpk, idex ) ! now i-velocity
	1105	CALL histwrite( id_i, "vomecrty", kt, vn , jpijpjjpk, idex ) ! now j-velocity
	1106	CALL histwrite( id_i, "vovecrtz", kt, wn , jpijpjjpk, idex ) ! now k-velocity
	1107	CALL histwrite( id_i, "sowaflup", kt, (emp-rnf ) , jpi*jpj , idex ) ! freshwater budget
	1108	CALL histwrite( id_i, "sohefldo", kt, qsr + qns , jpi*jpj , idex ) ! total heat flux
	1109	CALL histwrite( id_i, "soshfldo", kt, qsr , jpi*jpj , idex ) ! solar heat flux
	1110	CALL histwrite( id_i, "soicecov", kt, fr_i , jpi*jpj , idex ) ! ice fraction
	1111	CALL histwrite( id_i, "sozotaux", kt, utau , jpi*jpj , idex ) ! i-wind stress
	1112	CALL histwrite( id_i, "sometauy", kt, vtau , jpi*jpj , idex ) ! j-wind stress
[6204]	1113	IF( lk_vvl ) THEN
	1114	CALL histwrite( id_i, "vovvldep", kt, fsdept_n(:,:,:), jpijpjjpk, idex )! T-cell depth
	1115	CALL histwrite( id_i, "vovvle3t", kt, fse3t_n (:,:,:), jpijpjjpk, idex )! T-cell thickness
	1116	END IF
[3]	1117
	1118	! 3. Close the file
	1119	! -----------------
	1120	CALL histclo( id_i )
[1567]	1121	#if ! defined key_iomput && ! defined key_dimgout
[1561]	1122	IF( ninist /= 1 ) THEN
	1123	CALL histclo( nid_T )
	1124	CALL histclo( nid_U )
	1125	CALL histclo( nid_V )
	1126	CALL histclo( nid_W )
	1127	ENDIF
	1128	#endif
[3294]	1129
[3625]	1130	! IF( nn_timing == 1 ) CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
[3294]	1131	!
[3]	1132
	1133	END SUBROUTINE dia_wri_state
	1134	!!======================================================================
	1135	END MODULE diawri

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

Download in other formats: