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 NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90 @ 13151

Last change on this file since 13151 was 13151, checked in by gm, 4 years ago
result from merge with qco r12983
Property svn:keywords set to `Id`
File size: 54.2 KB

Rev	Line
[3]	1	MODULE diawri
	2	!!======================================================================
	3	!! * MODULE diawri *
	4	!! Ocean diagnostics : write ocean output files
	5	!!=====================================================================
[2528]	6	!! History : OPA ! 1991-03 (M.-A. Foujols) Original code
	7	!! 4.0 ! 1991-11 (G. Madec)
	8	!! ! 1992-06 (M. Imbard) correction restart file
	9	!! ! 1992-07 (M. Imbard) split into diawri and rstwri
	10	!! ! 1993-03 (M. Imbard) suppress writibm
	11	!! ! 1998-01 (C. Levy) NETCDF format using ioipsl INTERFACE
	12	!! ! 1999-02 (E. Guilyardi) name of netCDF files + variables
	13	!! 8.2 ! 2000-06 (M. Imbard) Original code (diabort.F)
	14	!! NEMO 1.0 ! 2002-06 (A.Bozec, E. Durand) Original code (diainit.F)
	15	!! - ! 2002-09 (G. Madec) F90: Free form and module
	16	!! - ! 2002-12 (G. Madec) merge of diabort and diainit, F90
	17	!! ! 2005-11 (V. Garnier) Surface pressure gradient organization
	18	!! 3.2 ! 2008-11 (B. Lemaire) creation from old diawri
[5836]	19	!! 3.7 ! 2014-01 (G. Madec) remove eddy induced velocity from no-IOM output
	20	!! ! change name of output variables in dia_wri_state
[2528]	21	!!----------------------------------------------------------------------
[3]	22
	23	!!----------------------------------------------------------------------
[2528]	24	!! dia_wri : create the standart output files
	25	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
	26	!!----------------------------------------------------------------------
[9019]	27	USE oce ! ocean dynamics and tracers
[12377]	28	USE isf_oce
	29	USE isfcpl
	30	USE abl ! abl variables in case ln_abl = .true.
[9019]	31	USE dom_oce ! ocean space and time domain
	32	USE phycst ! physical constants
	33	USE dianam ! build name of file (routine)
	34	USE diahth ! thermocline diagnostics
	35	USE dynadv , ONLY: ln_dynadv_vec
	36	USE icb_oce ! Icebergs
	37	USE icbdia ! Iceberg budgets
	38	USE ldftra ! lateral physics: eddy diffusivity coef.
	39	USE ldfdyn ! lateral physics: eddy viscosity coef.
	40	USE sbc_oce ! Surface boundary condition: ocean fields
	41	USE sbc_ice ! Surface boundary condition: ice fields
	42	USE sbcssr ! restoring term toward SST/SSS climatology
	43	USE sbcwave ! wave parameters
	44	USE wet_dry ! wetting and drying
	45	USE zdf_oce ! ocean vertical physics
	46	USE zdfdrg ! ocean vertical physics: top/bottom friction
	47	USE zdfmxl ! mixed layer
[6140]	48	!
[9019]	49	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	50	USE in_out_manager ! I/O manager
	51	USE dia25h ! 25h Mean output
	52	USE iom !
	53	USE ioipsl !
[5463]	54
[9570]	55	#if defined key_si3
[10425]	56	USE ice
[9019]	57	USE icewri
[1482]	58	#endif
[2715]	59	USE lib_mpp ! MPP library
[3294]	60	USE timing ! preformance summary
[12377]	61	USE diu_bulk ! diurnal warm layer
	62	USE diu_coolskin ! Cool skin
[2528]	63
[3]	64	IMPLICIT NONE
	65	PRIVATE
	66
[2528]	67	PUBLIC dia_wri ! routines called by step.F90
	68	PUBLIC dia_wri_state
[2715]	69	PUBLIC dia_wri_alloc ! Called by nemogcm module
[12377]	70	#if ! defined key_iomput
	71	PUBLIC dia_wri_alloc_abl ! Called by sbcabl module (if ln_abl = .true.)
	72	#endif
[2528]	73	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[3609]	74	INTEGER :: nb_T , ndim_bT ! grid_T file
[2528]	75	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	76	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	77	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
[12377]	78	INTEGER :: nid_A, nz_A, nh_A, ndim_A, ndim_hA ! grid_ABL file
[2528]	79	INTEGER :: ndex(1) ! ???
[2715]	80	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
[12377]	81	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL
[2715]	82	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]	83	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	84
	85	!! * Substitutions
[12377]	86	# include "do_loop_substitute.h90"
[13151]	87	!!st12
	88	# include "domzgr_substitute.h90"
[3]	89	!!----------------------------------------------------------------------
[9598]	90	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	91	!! $Id$
[10068]	92	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	93	!!----------------------------------------------------------------------
	94	CONTAINS
	95
[6140]	96	#if defined key_iomput
[3]	97	!!----------------------------------------------------------------------
[2528]	98	!! 'key_iomput' use IOM library
	99	!!----------------------------------------------------------------------
[9652]	100	INTEGER FUNCTION dia_wri_alloc()
	101	!
	102	dia_wri_alloc = 0
	103	!
	104	END FUNCTION dia_wri_alloc
[2715]	105
[9652]	106
[12377]	107	SUBROUTINE dia_wri( kt, Kmm )
[1482]	108	!!---------------------------------------------------------------------
	109	!! * ROUTINE dia_wri *
	110	!!
	111	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	112	!! NETCDF format is used by default
	113	!!
	114	!! ** Method : use iom_put
	115	!!----------------------------------------------------------------------
	116	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[12377]	117	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
[1756]	118	!!
[9019]	119	INTEGER :: ji, jj, jk ! dummy loop indices
	120	INTEGER :: ikbot ! local integer
	121	REAL(wp):: zztmp , zztmpx ! local scalar
	122	REAL(wp):: zztmp2, zztmpy ! - -
[12614]	123	REAL(wp):: ze3 ! - -
[9019]	124	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
	125	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
[1482]	126	!!----------------------------------------------------------------------
	127	!
[9124]	128	IF( ln_timing ) CALL timing_start('dia_wri')
[3294]	129	!
[1482]	130	! Output the initial state and forcings
	131	IF( ninist == 1 ) THEN
[12377]	132	CALL dia_wri_state( Kmm, 'output.init' )
[1482]	133	ninist = 0
	134	ENDIF
[3]	135
[6351]	136	! Output of initial vertical scale factor
	137	CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10114]	138	CALL iom_put("e3u_0", e3u_0(:,:,:) )
	139	CALL iom_put("e3v_0", e3v_0(:,:,:) )
[6351]	140	!
[13151]	141	!!st13
	142	#if ! defined key_qco
	143	IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN ! time-varying e3t
	144	DO jk = 1, jpk
	145	z3d(:,:,jk) = e3t(:,:,jk,Kmm)
	146	END DO
	147	CALL iom_put( "e3t" , z3d(:,:,:) )
	148	CALL iom_put( "e3tdef" , ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
	149	ENDIF
	150	IF ( iom_use("e3u") ) THEN ! time-varying e3u
	151	DO jk = 1, jpk
	152	z3d(:,:,jk) = e3u(:,:,jk,Kmm)
	153	END DO
	154	CALL iom_put( "e3u" , z3d(:,:,:) )
	155	ENDIF
	156	IF ( iom_use("e3v") ) THEN ! time-varying e3v
	157	DO jk = 1, jpk
	158	z3d(:,:,jk) = e3v(:,:,jk,Kmm)
	159	END DO
	160	CALL iom_put( "e3v" , z3d(:,:,:) )
	161	ENDIF
	162	IF ( iom_use("e3w") ) THEN ! time-varying e3w
	163	DO jk = 1, jpk
	164	z3d(:,:,jk) = e3w(:,:,jk,Kmm)
	165	END DO
	166	CALL iom_put( "e3w" , z3d(:,:,:) )
	167	ENDIF
	168	#endif
	169	!!st13
[9023]	170	IF( ll_wd ) THEN
[12377]	171	CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)
[9023]	172	ELSE
[12614]	173	CALL iom_put( "ssh" , ssh(:,:,Kmm) ) ! sea surface height
[9023]	174	ENDIF
	175
[12614]	176	!!an
	177	IF( iom_use("ht") ) & ! water column at t-point
	178	CALL iom_put( "ht" , ht_0(:,:) + ssh(:,:,Kmm) )
	179	!
	180	IF( iom_use("hu") ) THEN ! water column at u-point
	181	z2d(:,:) = 0._wp
	182	DO_2D_10_10
	183	z2d(ji,jj) = 0.5_wp * ( e3t(ji ,jj,1,Kmm) * e1e2t(ji ,jj) &
	184	& + e3t(ji+1,jj,1,Kmm) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj)
[12377]	185	END_2D
[12614]	186	CALL lbc_lnk( 'diawri', z2d, 'U', 1._wp )
	187	CALL iom_put( "hu", z2d )
[5107]	188	ENDIF
[12614]	189	!
	190	IF( iom_use("hv") ) THEN ! water column at v-point
	191	z2d(:,:) = 0._wp
	192	DO_2D_10_10
	193	z2d(ji,jj) = 0.5_wp * ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) &
	194	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) ) * r1_e1e2v(ji,jj)
[12377]	195	END_2D
[12614]	196	CALL lbc_lnk( 'diawri', z2d, 'V', 1._wp )
	197	CALL iom_put( "hv", z2d )
	198	ENDIF
	199	!
	200	IF( iom_use("hf") ) THEN ! water column at f-point
	201	z2d(:,:) = 0._wp
	202	DO_2D_10_10
	203	z2d(ji,jj) = 0.25_wp * ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	204	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	205	END_2D
[12667]	206	z2d(:,:) = z2d(:,:) * ssfmask(:,:)
[12614]	207	CALL lbc_lnk( 'diawri', z2d, 'F', 1._wp )
	208	CALL iom_put( "hf", z2d )
	209	ENDIF
	210	!!an
[5463]	211
[12614]	212	IF( iom_use("wetdep") ) & ! wet depth
	213	CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) )
	214
[5463]	215	IF ( iom_use("taubot") ) THEN ! bottom stress
[12489]	216	zztmp = rho0 * 0.25
[7753]	217	z2d(:,:) = 0._wp
[12377]	218	DO_2D_00_00
	219	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 &
	220	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 &
	221	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 &
	222	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2
	223	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
	224	!
	225	END_2D
[10425]	226	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[5463]	227	CALL iom_put( "taubot", z2d )
	228	ENDIF
[12614]	229	!
	230	CALL iom_put( "ssu" , uu(:,:,1,Kmm) ) ! surface i-current
	231	CALL iom_put( "ssv" , vv(:,:,1,Kmm) ) ! surface j-current
	232	CALL iom_put( "woce", ww ) ! vertical velocity
	233	!
[5107]	234
[12614]	235	IF ( iom_use("sKE") ) THEN ! surface kinetic energy at T point
	236	z2d(:,:) = 0._wp
	237	DO_2D_00_00
	238	z2d(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
	239	& + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm) &
	240	& + vv(ji,jj ,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,1,Kmm) &
	241	& + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm) ) &
	242	& * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * tmask(ji,jj,1)
[12377]	243	END_2D
[12614]	244	!
	245	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
	246	IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d )
	247
[4990]	248	ENDIF
[12614]	249
	250	!!an sKEf
	251	IF ( iom_use("sKEf") ) THEN ! surface kinetic energy at F point
	252	z2d(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
[12377]	253	DO_2D_00_00
[12614]	254	z2d(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
	255	& + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm) &
	256	& + vv(ji ,jj,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji ,jj) * e3v(ji ,jj,1,Kmm) &
	257	& + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm) ) &
	258	& * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
[12377]	259	END_2D
[12614]	260	!
	261	CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
	262	CALL iom_put( "sKEf", z2d )
[4990]	263	ENDIF
[12614]	264	!!an
[4840]	265	!
[12377]	266	CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence
[6351]	267	!
[12614]	268	! Output of vorticity terms
	269	IF ( iom_use("relvor") .OR. iom_use("plavor") .OR. &
	270	& iom_use("relpotvor") .OR. iom_use("abspotvor") .OR. &
	271	& iom_use("Ens") ) THEN
	272	!
[9019]	273	z2d(:,:) = 0._wp
[12614]	274	ze3 = 0._wp
	275	DO_2D_10_10
	276	z2d(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
	277	& - e1u(ji ,jj+1) * uu(ji ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm) ) * r1_e1e2f(ji,jj)
	278	END_2D
	279	CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
	280	CALL iom_put( "relvor", z2d ) ! relative vorticity ( zeta )
	281	!
	282	CALL iom_put( "plavor", ff_f ) ! planetary vorticity ( f )
	283	!
	284	DO_2D_10_10
	285	ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	286	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	287	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	288	ELSE ; ze3 = 0._wp
	289	ENDIF
	290	z2d(ji,jj) = ze3 * z2d(ji,jj)
	291	END_2D
	292	CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
	293	CALL iom_put( "relpotvor", z2d ) ! relative potential vorticity (zeta/h)
	294	!
	295	DO_2D_10_10
	296	ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	297	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	298	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	299	ELSE ; ze3 = 0._wp
	300	ENDIF
	301	z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)
	302	END_2D
	303	CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
	304	CALL iom_put( "abspotvor", z2d ) ! absolute potential vorticity ( q )
	305	!
	306	DO_2D_10_10
	307	z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj)
	308	END_2D
	309	CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
	310	CALL iom_put( "Ens", z2d ) ! potential enstrophy ( 1/2*q2 )
	311	!
[4990]	312	ENDIF
	313
[2528]	314	!
[9124]	315	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	316	!
[1482]	317	END SUBROUTINE dia_wri
	318
	319	#else
[2528]	320	!!----------------------------------------------------------------------
	321	!! Default option use IOIPSL library
	322	!!----------------------------------------------------------------------
	323
[9652]	324	INTEGER FUNCTION dia_wri_alloc()
	325	!!----------------------------------------------------------------------
	326	INTEGER, DIMENSION(2) :: ierr
	327	!!----------------------------------------------------------------------
[12493]	328	IF( nn_write == -1 ) THEN
	329	dia_wri_alloc = 0
	330	ELSE
	331	ierr = 0
	332	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	333	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	334	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
[9652]	335	!
[12493]	336	dia_wri_alloc = MAXVAL(ierr)
	337	CALL mpp_sum( 'diawri', dia_wri_alloc )
	338	!
	339	ENDIF
[9652]	340	!
	341	END FUNCTION dia_wri_alloc
[12377]	342
	343	INTEGER FUNCTION dia_wri_alloc_abl()
	344	!!----------------------------------------------------------------------
	345	ALLOCATE( ndex_hA(jpijpj), ndex_A (jpijpj*jpkam1), STAT=dia_wri_alloc_abl)
	346	CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
	347	!
	348	END FUNCTION dia_wri_alloc_abl
[9652]	349
	350
[12377]	351	SUBROUTINE dia_wri( kt, Kmm )
[3]	352	!!---------------------------------------------------------------------
	353	!! * ROUTINE dia_wri *
	354	!!
	355	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	356	!! NETCDF format is used by default
	357	!!
	358	!! ** Method : At the beginning of the first time step (nit000),
	359	!! define all the NETCDF files and fields
	360	!! At each time step call histdef to compute the mean if ncessary
[11536]	361	!! Each nn_write time step, output the instantaneous or mean fields
[3]	362	!!----------------------------------------------------------------------
[5836]	363	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[12377]	364	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
[5836]	365	!
[2528]	366	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	367	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	368	INTEGER :: inum = 11 ! temporary logical unit
[3294]	369	INTEGER :: ji, jj, jk ! dummy loop indices
	370	INTEGER :: ierr ! error code return from allocation
[2528]	371	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[12377]	372	INTEGER :: ipka ! ABL
[3609]	373	INTEGER :: jn, ierror ! local integers
[6140]	374	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	375	!
[9019]	376	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
[13151]	377	!!st14
	378	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept ! 3D workspace
[12377]	379	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace
[3]	380	!!----------------------------------------------------------------------
[1482]	381	!
[9019]	382	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[12377]	383	CALL dia_wri_state( Kmm, 'output.init' )
[1482]	384	ninist = 0
	385	ENDIF
	386	!
[11536]	387	IF( nn_write == -1 ) RETURN ! we will never do any output
	388	!
	389	IF( ln_timing ) CALL timing_start('dia_wri')
	390	!
[3]	391	! 0. Initialisation
	392	! -----------------
[632]	393
[9019]	394	ll_print = .FALSE. ! local variable for debugging
[3]	395	ll_print = ll_print .AND. lwp
	396
	397	! Define frequency of output and means
[5566]	398	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]	399	#if defined key_diainstant
[12489]	400	zsto = nn_write * rn_Dt
[1312]	401	clop = "inst("//TRIM(clop)//")"
[3]	402	#else
[12489]	403	zsto=rn_Dt
[1312]	404	clop = "ave("//TRIM(clop)//")"
[3]	405	#endif
[12489]	406	zout = nn_write * rn_Dt
	407	zmax = ( nitend - nit000 + 1 ) * rn_Dt
[3]	408
	409	! Define indices of the horizontal output zoom and vertical limit storage
	410	iimi = 1 ; iima = jpi
	411	ijmi = 1 ; ijma = jpj
	412	ipk = jpk
[12377]	413	IF(ln_abl) ipka = jpkam1
[3]	414
	415	! define time axis
[1334]	416	it = kt
	417	itmod = kt - nit000 + 1
[13151]	418	!!st15
	419	! store e3t for subsitute
	420	DO jk = 1, jpk
	421	ze3t (:,:,jk) = e3t (:,:,jk,Kmm)
	422	zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
	423	END DO
	424	!!st15
[3]	425
	426	! 1. Define NETCDF files and fields at beginning of first time step
	427	! -----------------------------------------------------------------
	428
	429	IF( kt == nit000 ) THEN
	430
	431	! Define the NETCDF files (one per grid)
[632]	432
[3]	433	! Compute julian date from starting date of the run
[12489]	434	CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
[1309]	435	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	436	IF(lwp)WRITE(numout,*)
	437	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	438	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	439	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	440	' limit storage in depth = ', ipk
	441
	442	! WRITE root name in date.file for use by postpro
[1581]	443	IF(lwp) THEN
[11536]	444	CALL dia_nam( clhstnam, nn_write,' ' )
[1581]	445	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	446	WRITE(inum,*) clhstnam
	447	CLOSE(inum)
	448	ENDIF
[632]	449
[3]	450	! Define the T grid FILE ( nid_T )
[632]	451
[11536]	452	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[3]	453	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	454	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	455	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	456	& nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	457	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	458	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	459	! ! Index of ocean points
	460	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	461	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	462	!
	463	IF( ln_icebergs ) THEN
	464	!
	465	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	466	!! that routine is called from nemogcm, so do it here immediately before its needed
	467	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	468	CALL mpp_sum( 'diawri', ierror )
[3609]	469	IF( ierror /= 0 ) THEN
	470	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	471	RETURN
	472	ENDIF
	473	!
	474	!! iceberg vertical coordinate is class number
	475	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	476	& "number", nclasses, class_num, nb_T )
	477	!
	478	!! each class just needs the surface index pattern
	479	ndim_bT = 3
	480	DO jn = 1,nclasses
	481	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	482	ENDDO
	483	!
	484	ENDIF
[3]	485
	486	! Define the U grid FILE ( nid_U )
	487
[11536]	488	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[3]	489	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	490	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	491	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	492	& nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	493	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	494	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	495	! ! Index of ocean points
	496	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	497	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	498
	499	! Define the V grid FILE ( nid_V )
	500
[11536]	501	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[3]	502	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	503	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	504	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	505	& nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	506	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	507	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	508	! ! Index of ocean points
	509	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	510	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	511
	512	! Define the W grid FILE ( nid_W )
	513
[11536]	514	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[3]	515	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	516	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	517	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	518	& nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	519	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	520	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	521
[12377]	522	IF( ln_abl ) THEN
	523	! Define the ABL grid FILE ( nid_A )
	524	CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
	525	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	526	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	527	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	528	& nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
[12377]	529	CALL histvert( nid_A, "ght_abl", "Vertical T levels", & ! Vertical grid: gdept
	530	& "m", ipka, ght_abl(2:jpka), nz_A, "up" )
	531	! ! Index of ocean points
	532	ALLOCATE( zw3d_abl(jpi,jpj,ipka) )
	533	zw3d_abl(:,:,:) = 1._wp
	534	CALL wheneq( jpijpjipka, zw3d_abl, 1, 1., ndex_A , ndim_A ) ! volume
	535	CALL wheneq( jpi*jpj , zw3d_abl, 1, 1., ndex_hA, ndim_hA ) ! surface
	536	DEALLOCATE(zw3d_abl)
	537	ENDIF
[632]	538
[3]	539	! Declare all the output fields as NETCDF variables
	540
	541	! !!! nid_T : 3D
	542	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	543	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	544	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	545	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]	546	IF( .NOT.ln_linssh ) THEN
[13151]	547	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! ze3t(:,:,:,Kmm)
[4292]	548	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[13151]	549	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! ze3t(:,:,:,Kmm)
[4292]	550	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[13151]	551	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! ze3t(:,:,:,Kmm)
[4292]	552	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	553	ENDIF
[3]	554	! !!! nid_T : 2D
	555	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	556	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	557	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	558	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	559	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	560	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	561	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	562	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	563	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	564	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	565	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	566	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[6140]	567	IF( ln_linssh ) THEN
[12377]	568	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * ts(:,:,1,jp_tem,Kmm)
[3625]	569	& , "KgC/m2/s", & ! sosst_cd
[4292]	570	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[12377]	571	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * ts(:,:,1,jp_sal,Kmm)
[3625]	572	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	573	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	574	ENDIF
[888]	575	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	576	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	577	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	578	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	579	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	580	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	581	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	582	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	583	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	584	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	585	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	586	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[12377]	587	!
	588	IF( ln_abl ) THEN
	589	CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl
	590	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	591	CALL histdef( nid_A, "q_abl", "Humidity" , "kg/kg" , & ! q_abl
	592	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	593	CALL histdef( nid_A, "u_abl", "Atmospheric U-wind " , "m/s" , & ! u_abl
	594	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	595	CALL histdef( nid_A, "v_abl", "Atmospheric V-wind " , "m/s" , & ! v_abl
	596	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	597	CALL histdef( nid_A, "tke_abl", "Atmospheric TKE " , "m2/s2" , & ! tke_abl
	598	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	599	CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s" , & ! avm_abl
	600	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	601	CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2", & ! avt_abl
	602	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	603	CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height " , "m", & ! pblh
	604	& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	605	#if defined key_si3
	606	CALL histdef( nid_A, "oce_frac", "Fraction of open ocean" , " ", & ! ato_i
	607	& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	608	#endif
	609	CALL histend( nid_A, snc4chunks=snc4set )
	610	ENDIF
	611	!
[3609]	612	IF( ln_icebergs ) THEN
	613	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	614	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	615	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	616	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	617	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	618	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	619	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	620	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	621	IF( ln_bergdia ) THEN
	622	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	623	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	624	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	625	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	626	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	627	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	628	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	629	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	630	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	631	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	632	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	633	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	634	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	635	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	636	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	637	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	638	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	639	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	640	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	641	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	642	ENDIF
	643	ENDIF
	644
[11536]	645	IF( ln_ssr ) THEN
[4990]	646	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	647	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	648	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	649	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	650	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	651	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	652	ENDIF
[11536]	653
[3]	654	clmx ="l_max(only(x))" ! max index on a period
[5836]	655	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	656	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
[3]	657	#if defined key_diahth
	658	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	659	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	660	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	661	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	662	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	663	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[7646]	664	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
[3]	665	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	666	#endif
	667
[2528]	668	CALL histend( nid_T, snc4chunks=snc4set )
[3]	669
	670	! !!! nid_U : 3D
[12377]	671	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! uu(:,:,:,Kmm)
[3]	672	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]	673	IF( ln_wave .AND. ln_sdw) THEN
	674	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	675	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	676	ENDIF
[3]	677	! !!! nid_U : 2D
[888]	678	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	679	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	680
[2528]	681	CALL histend( nid_U, snc4chunks=snc4set )
[3]	682
	683	! !!! nid_V : 3D
[12377]	684	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vv(:,:,:,Kmm)
[3]	685	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]	686	IF( ln_wave .AND. ln_sdw) THEN
	687	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	688	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	689	ENDIF
[3]	690	! !!! nid_V : 2D
[888]	691	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	692	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	693
[2528]	694	CALL histend( nid_V, snc4chunks=snc4set )
[3]	695
	696	! !!! nid_W : 3D
[12377]	697	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! ww
[3]	698	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	699	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	700	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]	701	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
[255]	702	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	703
[9019]	704	IF( ln_zdfddm ) THEN
[3]	705	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	706	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	707	ENDIF
[7646]	708
	709	IF( ln_wave .AND. ln_sdw) THEN
	710	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	711	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	712	ENDIF
[3]	713	! !!! nid_W : 2D
[2528]	714	CALL histend( nid_W, snc4chunks=snc4set )
[3]	715
	716	IF(lwp) WRITE(numout,*)
	717	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	718	IF(ll_print) CALL FLUSH(numout )
	719
	720	ENDIF
	721
	722	! 2. Start writing data
	723	! ---------------------
	724
[4292]	725	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	726	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	727	! donne le nombre d'elements, et ndex la liste des indices a sortir
	728
[11536]	729	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
[3]	730	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	731	WRITE(numout,*) '~~~~~~ '
	732	ENDIF
[13151]	733	!!st16
[6140]	734	IF( .NOT.ln_linssh ) THEN
[13151]	735	CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T ) ! heat content
	736	CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T ) ! salt content
	737	CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	738	CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
	739	!!st16
[4292]	740	ELSE
[12377]	741	CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature
	742	CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T ) ! salinity
	743	CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT ) ! sea surface temperature
	744	CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT ) ! sea surface salinity
[4292]	745	ENDIF
[6140]	746	IF( .NOT.ln_linssh ) THEN
[13151]	747	!!st17 if ! defined key_qco
	748	zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
	749	CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:) , ndim_T , ndex_T ) ! level thickness
	750	CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T ) ! t-point depth
	751	!!st17
[4292]	752	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	753	ENDIF
[12377]	754	CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm) , ndim_hT, ndex_hT ) ! sea surface height
[2528]	755	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	756	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	757	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	758	! (includes virtual salt flux beneath ice
	759	! in linear free surface case)
[6140]	760	IF( ln_linssh ) THEN
[12377]	761	zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
[4292]	762	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
[12377]	763	zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
[4292]	764	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	765	ENDIF
[888]	766	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	767	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	768	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	769	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	770	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	771	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[12377]	772	!
	773	IF( ln_abl ) THEN
	774	ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
	775	IF( ln_mskland ) THEN
	776	DO jk=1,jpka
	777	zw3d_abl(:,:,jk) = tmask(:,:,1)
	778	END DO
	779	ELSE
	780	zw3d_abl(:,:,:) = 1._wp
	781	ENDIF
	782	CALL histwrite( nid_A, "pblh" , it, pblh(:,:) *zw3d_abl(:,:,1 ), ndim_hA, ndex_hA ) ! pblh
	783	CALL histwrite( nid_A, "u_abl" , it, u_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! u_abl
	784	CALL histwrite( nid_A, "v_abl" , it, v_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! v_abl
	785	CALL histwrite( nid_A, "t_abl" , it, tq_abl (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! t_abl
	786	CALL histwrite( nid_A, "q_abl" , it, tq_abl (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! q_abl
	787	CALL histwrite( nid_A, "tke_abl", it, tke_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! tke_abl
	788	CALL histwrite( nid_A, "avm_abl", it, avm_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avm_abl
	789	CALL histwrite( nid_A, "avt_abl", it, avt_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avt_abl
	790	#if defined key_si3
	791	CALL histwrite( nid_A, "oce_frac" , it, ato_i(:,:) , ndim_hA, ndex_hA ) ! ato_i
	792	#endif
	793	DEALLOCATE(zw3d_abl)
	794	ENDIF
	795	!
[3609]	796	IF( ln_icebergs ) THEN
	797	!
	798	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	799	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	800	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	801	!
	802	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	803	!
	804	IF( ln_bergdia ) THEN
	805	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	806	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	807	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	808	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	809	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	810	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	811	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	812	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	813	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	814	!
	815	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	816	ENDIF
	817	ENDIF
	818
[11536]	819	IF( ln_ssr ) THEN
[4990]	820	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	821	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[12377]	822	zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
[4990]	823	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	824	ENDIF
	825	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	826	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	827
	828	#if defined key_diahth
	829	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	830	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	831	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	832	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	833	#endif
[888]	834
[12377]	835	CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U ) ! i-current
[888]	836	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	837
[12377]	838	CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V ) ! j-current
[888]	839	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	840
[11418]	841	IF( ln_zad_Aimp ) THEN
[12377]	842	CALL histwrite( nid_W, "vovecrtz", it, ww + wi , ndim_T, ndex_T ) ! vert. current
[11418]	843	ELSE
[12377]	844	CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current
[11418]	845	ENDIF
[3]	846	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	847	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	848	IF( ln_zdfddm ) THEN
	849	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	850	ENDIF
	851
[7646]	852	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	853	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	854	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	855	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	856	ENDIF
	857
[1318]	858	! 3. Close all files
	859	! ---------------------------------------
[1561]	860	IF( kt == nitend ) THEN
[3]	861	CALL histclo( nid_T )
	862	CALL histclo( nid_U )
	863	CALL histclo( nid_V )
	864	CALL histclo( nid_W )
[12377]	865	IF(ln_abl) CALL histclo( nid_A )
[3]	866	ENDIF
[2528]	867	!
[9124]	868	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	869	!
[3]	870	END SUBROUTINE dia_wri
[1567]	871	#endif
	872
[12377]	873	SUBROUTINE dia_wri_state( Kmm, cdfile_name )
[3]	874	!!---------------------------------------------------------------------
	875	!! * ROUTINE dia_wri_state *
	876	!!
	877	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	878	!! the instantaneous ocean state and forcing fields.
	879	!! Used to find errors in the initial state or save the last
	880	!! ocean state in case of abnormal end of a simulation
	881	!!
	882	!! ** Method : NetCDF files using ioipsl
	883	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	884	!! File 'output.abort.nc' is created in case of abnormal job end
	885	!!----------------------------------------------------------------------
[12377]	886	INTEGER , INTENT( in ) :: Kmm ! time level index
[1334]	887	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	888	!!
[12377]	889	INTEGER :: inum, jk
[13151]	890	!!st18 TBR
	891	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept ! 3D workspace !!st patch to use substitution
[3]	892	!!----------------------------------------------------------------------
[3294]	893	!
[648]	894	IF(lwp) WRITE(numout,*)
	895	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	896	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
[10425]	897	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
[13151]	898	!!st19 TBR
	899	DO jk = 1, jpk
	900	ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
	901	zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
	902	END DO
	903	!!st19
[9570]	904	#if defined key_si3
[10425]	905	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
[1482]	906	#else
[10425]	907	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
[1482]	908	#endif
[3]	909
[12377]	910	CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) ) ! now temperature
	911	CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) ) ! now salinity
	912	CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm) ) ! sea surface height
	913	CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm) ) ! now i-velocity
	914	CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm) ) ! now j-velocity
[11418]	915	IF( ln_zad_Aimp ) THEN
[12377]	916	CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi ) ! now k-velocity
[11418]	917	ELSE
[12377]	918	CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity
[11418]	919	ENDIF
[12377]	920	CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep ) ! now k-velocity
[13151]	921	CALL iom_rstput( 0, 0, inum, 'ht' , ht(:,:) ) ! now water column height
[12377]	922
	923	IF ( ln_isf ) THEN
	924	IF (ln_isfcav_mlt) THEN
	925	CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav ) ! now k-velocity
	926	CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav ) ! now k-velocity
	927	CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav ) ! now k-velocity
[13151]	928	CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) ) ! now k-velocity
	929	CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) ) ! now k-velocity
	930	CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 )
[12377]	931	END IF
	932	IF (ln_isfpar_mlt) THEN
[13151]	933	CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) ) ! now k-velocity
[12377]	934	CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par ) ! now k-velocity
	935	CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par ) ! now k-velocity
	936	CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par ) ! now k-velocity
[13151]	937	CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) ) ! now k-velocity
	938	CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) ) ! now k-velocity
	939	CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 )
[12377]	940	END IF
	941	END IF
[13151]	942	!
[7646]	943	IF( ALLOCATED(ahtu) ) THEN
[10425]	944	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	945	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	946	ENDIF
	947	IF( ALLOCATED(ahmt) ) THEN
[10425]	948	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	949	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	950	ENDIF
[10425]	951	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	952	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	953	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	954	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	955	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	956	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[13151]	957	!!st20 TBR
	958	IF( .NOT.ln_linssh ) THEN
	959	CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept ) ! T-cell depth
	960	CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t ) ! T-cell thickness
[10425]	961	END IF
[7646]	962	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	963	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	964	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	965	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	966	ENDIF
[12377]	967	IF ( ln_abl ) THEN
	968	CALL iom_rstput ( 0, 0, inum, "uz1_abl", u_abl(:,:,2,nt_a ) ) ! now first level i-wind
	969	CALL iom_rstput ( 0, 0, inum, "vz1_abl", v_abl(:,:,2,nt_a ) ) ! now first level j-wind
	970	CALL iom_rstput ( 0, 0, inum, "tz1_abl", tq_abl(:,:,2,nt_a,1) ) ! now first level temperature
	971	CALL iom_rstput ( 0, 0, inum, "qz1_abl", tq_abl(:,:,2,nt_a,2) ) ! now first level humidity
	972	ENDIF
[10425]	973
	974	#if defined key_si3
	975	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
	976	CALL ice_wri_state( inum )
[1561]	977	ENDIF
	978	#endif
[10425]	979	!
	980	CALL iom_close( inum )
[3294]	981	!
[3]	982	END SUBROUTINE dia_wri_state
[6140]	983
[3]	984	!!======================================================================
	985	END MODULE diawri

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

Download in other formats: