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/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DIA/diawri.F90 @ 11501

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

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

Download in other formats: