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

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

diawri.F90 in branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5282

Last change on this file since 5282 was 5282, checked in by diovino, 9 years ago
Dev. branch CMCC4_simplification ticket #1456
Property svn:keywords set to `Id`
File size: 55.5 KB

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

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

Download in other formats: