Context Navigation

diawri.F90 @ 15343

Last change on this file since 15343 was 15333, checked in by hadjt, 3 years ago

Adding Regional Means, but without XIOS or MLD.

Search on
!JT MLD
!JT IOM

in IOM/iom.F90 and DIA/diaregmean.F90

to see the code commented out.

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diawri.F90 @ 15343

Download in other formats: