Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

diawri.F90 @ 12494

Last change on this file since 12494 was 12494, checked in by smasson, 4 years ago
release-4.0-HEAD: minor optimisation in diawri, see #2393
Property svn:keywords set to `Id`
File size: 49.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
	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	!!----------------------------------------------------------------------
[12494]	408	IF( nn_write == -1 ) THEN
	409	dia_wri_alloc = 0
	410	ELSE
	411	ierr = 0
	412	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	413	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	414	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
[9652]	415	!
[12494]	416	dia_wri_alloc = MAXVAL(ierr)
	417	CALL mpp_sum( 'diawri', dia_wri_alloc )
	418	!
	419	ENDIF
[9652]	420	!
	421	END FUNCTION dia_wri_alloc
	422
	423
[1561]	424	SUBROUTINE dia_wri( kt )
[3]	425	!!---------------------------------------------------------------------
	426	!! * ROUTINE dia_wri *
	427	!!
	428	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	429	!! NETCDF format is used by default
	430	!!
	431	!! ** Method : At the beginning of the first time step (nit000),
	432	!! define all the NETCDF files and fields
	433	!! At each time step call histdef to compute the mean if ncessary
[11536]	434	!! Each nn_write time step, output the instantaneous or mean fields
[3]	435	!!----------------------------------------------------------------------
[5836]	436	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	437	!
[2528]	438	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	439	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	440	INTEGER :: inum = 11 ! temporary logical unit
[3294]	441	INTEGER :: ji, jj, jk ! dummy loop indices
	442	INTEGER :: ierr ! error code return from allocation
[2528]	443	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[3609]	444	INTEGER :: jn, ierror ! local integers
[6140]	445	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	446	!
[9019]	447	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
	448	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
[3]	449	!!----------------------------------------------------------------------
[1482]	450	!
[9019]	451	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[10425]	452	CALL dia_wri_state( 'output.init' )
[1482]	453	ninist = 0
	454	ENDIF
	455	!
[11536]	456	IF( nn_write == -1 ) RETURN ! we will never do any output
	457	!
	458	IF( ln_timing ) CALL timing_start('dia_wri')
	459	!
[3]	460	! 0. Initialisation
	461	! -----------------
[632]	462
[9019]	463	ll_print = .FALSE. ! local variable for debugging
[3]	464	ll_print = ll_print .AND. lwp
	465
	466	! Define frequency of output and means
[5566]	467	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]	468	#if defined key_diainstant
[11536]	469	zsto = nn_write * rdt
[1312]	470	clop = "inst("//TRIM(clop)//")"
[3]	471	#else
[6140]	472	zsto=rdt
[1312]	473	clop = "ave("//TRIM(clop)//")"
[3]	474	#endif
[11536]	475	zout = nn_write * rdt
[6140]	476	zmax = ( nitend - nit000 + 1 ) * rdt
[3]	477
	478	! Define indices of the horizontal output zoom and vertical limit storage
	479	iimi = 1 ; iima = jpi
	480	ijmi = 1 ; ijma = jpj
	481	ipk = jpk
	482
	483	! define time axis
[1334]	484	it = kt
	485	itmod = kt - nit000 + 1
[3]	486
	487
	488	! 1. Define NETCDF files and fields at beginning of first time step
	489	! -----------------------------------------------------------------
	490
	491	IF( kt == nit000 ) THEN
	492
	493	! Define the NETCDF files (one per grid)
[632]	494
[3]	495	! Compute julian date from starting date of the run
[1309]	496	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	497	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	498	IF(lwp)WRITE(numout,*)
	499	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	500	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	501	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	502	' limit storage in depth = ', ipk
	503
	504	! WRITE root name in date.file for use by postpro
[1581]	505	IF(lwp) THEN
[11536]	506	CALL dia_nam( clhstnam, nn_write,' ' )
[1581]	507	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	508	WRITE(inum,*) clhstnam
	509	CLOSE(inum)
	510	ENDIF
[632]	511
[3]	512	! Define the T grid FILE ( nid_T )
[632]	513
[11536]	514	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[3]	515	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	516	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	517	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	518	& nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	519	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	520	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	521	! ! Index of ocean points
	522	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	523	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	524	!
	525	IF( ln_icebergs ) THEN
	526	!
	527	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	528	!! that routine is called from nemogcm, so do it here immediately before its needed
	529	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	530	CALL mpp_sum( 'diawri', ierror )
[3609]	531	IF( ierror /= 0 ) THEN
	532	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	533	RETURN
	534	ENDIF
	535	!
	536	!! iceberg vertical coordinate is class number
	537	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	538	& "number", nclasses, class_num, nb_T )
	539	!
	540	!! each class just needs the surface index pattern
	541	ndim_bT = 3
	542	DO jn = 1,nclasses
	543	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	544	ENDDO
	545	!
	546	ENDIF
[3]	547
	548	! Define the U grid FILE ( nid_U )
	549
[11536]	550	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[3]	551	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	552	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	553	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	554	& nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	555	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	556	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	557	! ! Index of ocean points
	558	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	559	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	560
	561	! Define the V grid FILE ( nid_V )
	562
[11536]	563	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[3]	564	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	565	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	566	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	567	& nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	568	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	569	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	570	! ! Index of ocean points
	571	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	572	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	573
	574	! Define the W grid FILE ( nid_W )
	575
[11536]	576	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[3]	577	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	578	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	579	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	580	& nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	581	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	582	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	583
[632]	584
[3]	585	! Declare all the output fields as NETCDF variables
	586
	587	! !!! nid_T : 3D
	588	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	589	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	590	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	591	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]	592	IF( .NOT.ln_linssh ) THEN
[4292]	593	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	594	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	595	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	596	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	597	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	598	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	599	ENDIF
[3]	600	! !!! nid_T : 2D
	601	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	602	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	603	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	604	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	605	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	606	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	607	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	608	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	609	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	610	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	611	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	612	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[6140]	613	IF( ln_linssh ) THEN
[4292]	614	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
[3625]	615	& , "KgC/m2/s", & ! sosst_cd
[4292]	616	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	617	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
[3625]	618	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	619	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	620	ENDIF
[888]	621	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	622	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	623	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	624	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	625	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	626	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	627	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	628	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	629	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	630	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	631	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	632	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3609]	633	!
	634	IF( ln_icebergs ) THEN
	635	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	636	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	637	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	638	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	639	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	640	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	641	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	642	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	643	IF( ln_bergdia ) THEN
	644	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	645	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	646	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	647	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	648	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	649	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	650	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	651	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	652	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	653	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	654	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	655	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	656	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	657	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	658	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	659	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	660	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	661	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	662	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	663	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	664	ENDIF
	665	ENDIF
	666
[11536]	667	IF( ln_ssr ) THEN
[4990]	668	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	669	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	670	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	671	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	672	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	673	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	674	ENDIF
[11536]	675
[3]	676	clmx ="l_max(only(x))" ! max index on a period
[5836]	677	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	678	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
[3]	679	#if defined key_diahth
	680	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	681	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	682	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	683	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	684	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	685	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[7646]	686	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
[3]	687	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	688	#endif
	689
[2528]	690	CALL histend( nid_T, snc4chunks=snc4set )
[3]	691
	692	! !!! nid_U : 3D
	693	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	694	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]	695	IF( ln_wave .AND. ln_sdw) THEN
	696	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	697	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	698	ENDIF
[3]	699	! !!! nid_U : 2D
[888]	700	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	701	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	702
[2528]	703	CALL histend( nid_U, snc4chunks=snc4set )
[3]	704
	705	! !!! nid_V : 3D
	706	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	707	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]	708	IF( ln_wave .AND. ln_sdw) THEN
	709	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	710	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	711	ENDIF
[3]	712	! !!! nid_V : 2D
[888]	713	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	714	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	715
[2528]	716	CALL histend( nid_V, snc4chunks=snc4set )
[3]	717
	718	! !!! nid_W : 3D
	719	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	720	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	721	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	722	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]	723	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
[255]	724	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	725
[9019]	726	IF( ln_zdfddm ) THEN
[3]	727	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	728	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	729	ENDIF
[7646]	730
	731	IF( ln_wave .AND. ln_sdw) THEN
	732	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	733	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	734	ENDIF
[3]	735	! !!! nid_W : 2D
[2528]	736	CALL histend( nid_W, snc4chunks=snc4set )
[3]	737
	738	IF(lwp) WRITE(numout,*)
	739	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	740	IF(ll_print) CALL FLUSH(numout )
	741
	742	ENDIF
	743
	744	! 2. Start writing data
	745	! ---------------------
	746
[4292]	747	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	748	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	749	! donne le nombre d'elements, et ndex la liste des indices a sortir
	750
[11536]	751	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
[3]	752	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	753	WRITE(numout,*) '~~~~~~ '
	754	ENDIF
	755
[6140]	756	IF( .NOT.ln_linssh ) THEN
	757	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	758	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	759	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	760	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	761	ELSE
	762	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	763	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	764	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	765	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	766	ENDIF
[6140]	767	IF( .NOT.ln_linssh ) THEN
[7753]	768	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
[6140]	769	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	770	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
[4292]	771	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	772	ENDIF
[3]	773	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
[2528]	774	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	775	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	776	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	777	! (includes virtual salt flux beneath ice
	778	! in linear free surface case)
[6140]	779	IF( ln_linssh ) THEN
[7753]	780	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
[4292]	781	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
[7753]	782	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
[4292]	783	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	784	ENDIF
[888]	785	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	786	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	787	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	788	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	789	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	790	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[3609]	791	!
	792	IF( ln_icebergs ) THEN
	793	!
	794	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	795	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	796	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	797	!
	798	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	799	!
	800	IF( ln_bergdia ) THEN
	801	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	802	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	803	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	804	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	805	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	806	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	807	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	808	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	809	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	810	!
	811	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	812	ENDIF
	813	ENDIF
	814
[11536]	815	IF( ln_ssr ) THEN
[4990]	816	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	817	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[11536]	818	zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	819	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	820	ENDIF
	821	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	822	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	823
	824	#if defined key_diahth
	825	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	826	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	827	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	828	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	829	#endif
[888]	830
[3]	831	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
[888]	832	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	833
	834	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
[888]	835	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	836
[11418]	837	IF( ln_zad_Aimp ) THEN
	838	CALL histwrite( nid_W, "vovecrtz", it, wn + wi , ndim_T, ndex_T ) ! vert. current
	839	ELSE
	840	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
	841	ENDIF
[3]	842	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	843	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	844	IF( ln_zdfddm ) THEN
	845	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	846	ENDIF
	847
[7646]	848	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	849	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	850	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	851	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	852	ENDIF
	853
[1318]	854	! 3. Close all files
	855	! ---------------------------------------
[1561]	856	IF( kt == nitend ) THEN
[3]	857	CALL histclo( nid_T )
	858	CALL histclo( nid_U )
	859	CALL histclo( nid_V )
	860	CALL histclo( nid_W )
	861	ENDIF
[2528]	862	!
[9124]	863	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	864	!
[3]	865	END SUBROUTINE dia_wri
[1567]	866	#endif
	867
[10425]	868	SUBROUTINE dia_wri_state( cdfile_name )
[3]	869	!!---------------------------------------------------------------------
	870	!! * ROUTINE dia_wri_state *
	871	!!
	872	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	873	!! the instantaneous ocean state and forcing fields.
	874	!! Used to find errors in the initial state or save the last
	875	!! ocean state in case of abnormal end of a simulation
	876	!!
	877	!! ** Method : NetCDF files using ioipsl
	878	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	879	!! File 'output.abort.nc' is created in case of abnormal job end
	880	!!----------------------------------------------------------------------
[1334]	881	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	882	!!
	883	INTEGER :: inum
[3]	884	!!----------------------------------------------------------------------
[3294]	885	!
[648]	886	IF(lwp) WRITE(numout,*)
	887	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	888	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
[10425]	889	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
[648]	890
[9570]	891	#if defined key_si3
[10425]	892	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
[1482]	893	#else
[10425]	894	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
[1482]	895	#endif
[3]	896
[10425]	897	CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) ) ! now temperature
	898	CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) ) ! now salinity
	899	CALL iom_rstput( 0, 0, inum, 'sossheig', sshn ) ! sea surface height
	900	CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity
	901	CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity
[11418]	902	IF( ln_zad_Aimp ) THEN
	903	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi ) ! now k-velocity
	904	ELSE
	905	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity
	906	ENDIF
[7646]	907	IF( ALLOCATED(ahtu) ) THEN
[10425]	908	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	909	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	910	ENDIF
	911	IF( ALLOCATED(ahmt) ) THEN
[10425]	912	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	913	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	914	ENDIF
[10425]	915	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	916	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	917	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	918	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	919	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	920	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[6140]	921	IF( .NOT.ln_linssh ) THEN
[10425]	922	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
	923	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
	924	END IF
[7646]	925	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	926	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	927	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	928	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	929	ENDIF
[10425]	930
	931	#if defined key_si3
	932	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
	933	CALL ice_wri_state( inum )
[1561]	934	ENDIF
	935	#endif
[10425]	936	!
	937	CALL iom_close( inum )
[3294]	938	!
[3]	939	END SUBROUTINE dia_wri_state
[6140]	940
[3]	941	!!======================================================================
	942	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/releases/release-4.0-HEAD/src/OCE/DIA/diawri.F90 @ 12494

Download in other formats: