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/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5433

Last change on this file since 5433 was 5260, checked in by deazer, 9 years ago
Merged branch with Trunk at revision 5253. Checked with SETTE, passes modified iodef.xml for AMM12 experiment
Property svn:keywords set to `Id`
File size: 54.5 KB

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

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

Download in other formats: