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/releases/r4.0/r4.0-HEAD/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DIA/diawri.F90 @ 14719

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

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

Download in other formats: