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/UKMO/NEMO_4.0.4_GO8_package/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/DIA/diawri.F90 @ 14078

Last change on this file since 14078 was 14078, checked in by cguiavarch, 3 years ago
UKMO/NEMO_4.0.4_GO8_package: copy over changes from NEMO_4.0.3_GO8_package branch @13647.
File size: 50.8 KB

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

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

Download in other formats: