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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diawri.F90 in NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DIA/diawri.F90 @ 10358

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

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

Download in other formats: