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

Context Navigation

source: branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5446

Last change on this file since 5446 was 5446, checked in by davestorkey, 9 years ago
Update UKMO/dev_r5021_nn_etau_revision branch to revision 5442 of the trunk.
File size: 54.6 KB

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

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

Download in other formats: