Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

diawri.F90 @ 13049

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/OCE/DIA/diawri.F90 @ 13049

Download in other formats: