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

Context Navigation

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4792

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

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

Download in other formats: