Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

diawri.F90 @ 12658

Last change on this file since 12658 was 12658, checked in by cguiavarch, 4 years ago
UKMO/NEMO_4.0.2_mirror : Remove SVN keywords.
File size: 49.0 KB

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

Download in other formats: