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

Context Navigation

source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5104

Last change on this file since 5104 was 5104, checked in by mathiot, 9 years ago
correction of minor bug with isf + add some extra test on ln_isfcav + restore ssu/ssv
Property svn:keywords set to `Id`
File size: 54.6 KB

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

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

Download in other formats: