Context Navigation

diawri.F90 @ 9495

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

Download in other formats: