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

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

diawri.F90 in branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4901

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

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

Download in other formats: