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/ENHANCE-02_ISF_nemo/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90 @ 11876

Last change on this file since 11876 was 11876, checked in by mathiot, 4 years ago
ENHANCE-02_ISF_nemo: add timing in the main isf routine + various bug fixes + cosmetic changes
Property svn:keywords set to `Id`
File size: 51.2 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
[11395]	28	USE isf
[11521]	29	USE isfcpl
[9019]	30	USE dom_oce ! ocean space and time domain
	31	USE phycst ! physical constants
	32	USE dianam ! build name of file (routine)
	33	USE diahth ! thermocline diagnostics
	34	USE dynadv , ONLY: ln_dynadv_vec
	35	USE icb_oce ! Icebergs
	36	USE icbdia ! Iceberg budgets
	37	USE ldftra ! lateral physics: eddy diffusivity coef.
	38	USE ldfdyn ! lateral physics: eddy viscosity coef.
	39	USE sbc_oce ! Surface boundary condition: ocean fields
	40	USE sbc_ice ! Surface boundary condition: ice fields
	41	USE sbcssr ! restoring term toward SST/SSS climatology
	42	USE sbcwave ! wave parameters
	43	USE wet_dry ! wetting and drying
	44	USE zdf_oce ! ocean vertical physics
	45	USE zdfdrg ! ocean vertical physics: top/bottom friction
	46	USE zdfmxl ! mixed layer
[6140]	47	!
[9019]	48	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	49	USE in_out_manager ! I/O manager
	50	USE diatmb ! Top,middle,bottom output
	51	USE dia25h ! 25h Mean output
	52	USE iom !
	53	USE ioipsl !
[5463]	54
[9570]	55	#if defined key_si3
[10425]	56	USE ice
[9019]	57	USE icewri
[1482]	58	#endif
[2715]	59	USE lib_mpp ! MPP library
[3294]	60	USE timing ! preformance summary
[6140]	61	USE diurnal_bulk ! diurnal warm layer
	62	USE cool_skin ! Cool skin
[2528]	63
[3]	64	IMPLICIT NONE
	65	PRIVATE
	66
[2528]	67	PUBLIC dia_wri ! routines called by step.F90
	68	PUBLIC dia_wri_state
[2715]	69	PUBLIC dia_wri_alloc ! Called by nemogcm module
[3]	70
[2528]	71	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[3609]	72	INTEGER :: nb_T , ndim_bT ! grid_T file
[2528]	73	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	74	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	75	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
	76	INTEGER :: ndex(1) ! ???
[2715]	77	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
	78	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]	79	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	80
	81	!! * Substitutions
[1756]	82	# include "vectopt_loop_substitute.h90"
[3]	83	!!----------------------------------------------------------------------
[9598]	84	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	85	!! $Id$
[10068]	86	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	87	!!----------------------------------------------------------------------
	88	CONTAINS
	89
[6140]	90	#if defined key_iomput
[3]	91	!!----------------------------------------------------------------------
[2528]	92	!! 'key_iomput' use IOM library
	93	!!----------------------------------------------------------------------
[9652]	94	INTEGER FUNCTION dia_wri_alloc()
	95	!
	96	dia_wri_alloc = 0
	97	!
	98	END FUNCTION dia_wri_alloc
[2715]	99
[9652]	100
[1561]	101	SUBROUTINE dia_wri( kt )
[1482]	102	!!---------------------------------------------------------------------
	103	!! * ROUTINE dia_wri *
	104	!!
	105	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	106	!! NETCDF format is used by default
	107	!!
	108	!! ** Method : use iom_put
	109	!!----------------------------------------------------------------------
	110	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1756]	111	!!
[9019]	112	INTEGER :: ji, jj, jk ! dummy loop indices
	113	INTEGER :: ikbot ! local integer
	114	REAL(wp):: zztmp , zztmpx ! local scalar
	115	REAL(wp):: zztmp2, zztmpy ! - -
	116	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
	117	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
[1482]	118	!!----------------------------------------------------------------------
	119	!
[9124]	120	IF( ln_timing ) CALL timing_start('dia_wri')
[3294]	121	!
[1482]	122	! Output the initial state and forcings
	123	IF( ninist == 1 ) THEN
[10425]	124	CALL dia_wri_state( 'output.init' )
[1482]	125	ninist = 0
	126	ENDIF
[3]	127
[6351]	128	! Output of initial vertical scale factor
	129	CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10114]	130	CALL iom_put("e3u_0", e3u_0(:,:,:) )
	131	CALL iom_put("e3v_0", e3v_0(:,:,:) )
[6351]	132	!
[6387]	133	CALL iom_put( "e3t" , e3t_n(:,:,:) )
	134	CALL iom_put( "e3u" , e3u_n(:,:,:) )
	135	CALL iom_put( "e3v" , e3v_n(:,:,:) )
	136	CALL iom_put( "e3w" , e3w_n(:,:,:) )
[6351]	137	IF( iom_use("e3tdef") ) &
[6387]	138	CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
[5461]	139
[9023]	140	IF( ll_wd ) THEN
	141	CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)
	142	ELSE
	143	CALL iom_put( "ssh" , sshn ) ! sea surface height
	144	ENDIF
	145
[7646]	146	IF( iom_use("wetdep") ) & ! wet depth
[9023]	147	CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
[5107]	148
	149	CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature
	150	CALL iom_put( "sst", tsn(:,:,1,jp_tem) ) ! surface temperature
	151	IF ( iom_use("sbt") ) THEN
[4990]	152	DO jj = 1, jpj
	153	DO ji = 1, jpi
[9019]	154	ikbot = mbkt(ji,jj)
	155	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
[4990]	156	END DO
[5107]	157	END DO
	158	CALL iom_put( "sbt", z2d ) ! bottom temperature
	159	ENDIF
	160
	161	CALL iom_put( "soce", tsn(:,:,:,jp_sal) ) ! 3D salinity
	162	CALL iom_put( "sss", tsn(:,:,1,jp_sal) ) ! surface salinity
	163	IF ( iom_use("sbs") ) THEN
[4990]	164	DO jj = 1, jpj
	165	DO ji = 1, jpi
[9019]	166	ikbot = mbkt(ji,jj)
	167	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
[4990]	168	END DO
[5107]	169	END DO
	170	CALL iom_put( "sbs", z2d ) ! bottom salinity
	171	ENDIF
[5463]	172
	173	IF ( iom_use("taubot") ) THEN ! bottom stress
[9019]	174	zztmp = rau0 * 0.25
[7753]	175	z2d(:,:) = 0._wp
[5463]	176	DO jj = 2, jpjm1
	177	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	178	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * un(ji ,jj,mbku(ji ,jj)) )**2 &
	179	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj)) )**2 &
	180	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vn(ji,jj ,mbkv(ji,jj )) )**2 &
	181	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1)) )**2
	182	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
[5463]	183	!
[9019]	184	END DO
	185	END DO
[10425]	186	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[5463]	187	CALL iom_put( "taubot", z2d )
	188	ENDIF
[5107]	189
[9019]	190	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
	191	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
[5107]	192	IF ( iom_use("sbu") ) THEN
[4990]	193	DO jj = 1, jpj
	194	DO ji = 1, jpi
[9019]	195	ikbot = mbku(ji,jj)
	196	z2d(ji,jj) = un(ji,jj,ikbot)
[4990]	197	END DO
[5107]	198	END DO
	199	CALL iom_put( "sbu", z2d ) ! bottom i-current
	200	ENDIF
	201
[9019]	202	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
	203	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
[5107]	204	IF ( iom_use("sbv") ) THEN
[4990]	205	DO jj = 1, jpj
	206	DO ji = 1, jpi
[9019]	207	ikbot = mbkv(ji,jj)
	208	z2d(ji,jj) = vn(ji,jj,ikbot)
[4990]	209	END DO
[5107]	210	END DO
	211	CALL iom_put( "sbv", z2d ) ! bottom j-current
[4990]	212	ENDIF
[1482]	213
[5461]	214	CALL iom_put( "woce", wn ) ! vertical velocity
	215	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
	216	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[7753]	217	z2d(:,:) = rau0 * e1e2t(:,:)
[5461]	218	DO jk = 1, jpk
[7753]	219	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
[5461]	220	END DO
	221	CALL iom_put( "w_masstr" , z3d )
	222	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
	223	ENDIF
	224
[9019]	225	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	226	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
	227	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
[5107]	228
[9019]	229	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
	230	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
[6351]	231
[5107]	232	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
[4990]	233	DO jj = 2, jpjm1 ! sst gradient
	234	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	235	zztmp = tsn(ji,jj,1,jp_tem)
	236	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)
	237	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]	238	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	239	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
	240	END DO
[1756]	241	END DO
[10425]	242	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[9019]	243	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
[7753]	244	z2d(:,:) = SQRT( z2d(:,:) )
[9019]	245	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
[4990]	246	ENDIF
	247
[9019]	248	! heat and salt contents
[4990]	249	IF( iom_use("heatc") ) THEN
[7753]	250	z2d(:,:) = 0._wp
[4990]	251	DO jk = 1, jpkm1
[5107]	252	DO jj = 1, jpj
	253	DO ji = 1, jpi
[6140]	254	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
[4990]	255	END DO
[4761]	256	END DO
	257	END DO
[9019]	258	CALL iom_put( "heatc", rau0_rcp * z2d ) ! vertically integrated heat content (J/m2)
[4990]	259	ENDIF
	260
	261	IF( iom_use("saltc") ) THEN
[7753]	262	z2d(:,:) = 0._wp
[4990]	263	DO jk = 1, jpkm1
[5107]	264	DO jj = 1, jpj
	265	DO ji = 1, jpi
[6140]	266	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
[4990]	267	END DO
	268	END DO
	269	END DO
[9019]	270	CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
[4990]	271	ENDIF
[4840]	272	!
[4990]	273	IF ( iom_use("eken") ) THEN
[9399]	274	z3d(:,:,jpk) = 0._wp
[4990]	275	DO jk = 1, jpkm1
	276	DO jj = 2, jpjm1
	277	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	278	zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
	279	z3d(ji,jj,jk) = zztmp * ( un(ji-1,jj,jk)*2 e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) &
	280	& + un(ji ,jj,jk)*2 e2u(ji ,jj) * e3u_n(ji ,jj,jk) &
	281	& + vn(ji,jj-1,jk)*2 e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) &
	282	& + vn(ji,jj ,jk)*2 e1v(ji,jj ) * e3v_n(ji,jj ,jk) )
	283	END DO
	284	END DO
	285	END DO
[10425]	286	CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
[9019]	287	CALL iom_put( "eken", z3d ) ! kinetic energy
[4990]	288	ENDIF
[6351]	289	!
	290	CALL iom_put( "hdiv", hdivn ) ! Horizontal divergence
	291	!
[7646]	292	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
[7753]	293	z3d(:,:,jpk) = 0.e0
	294	z2d(:,:) = 0.e0
[1756]	295	DO jk = 1, jpkm1
[7753]	296	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
	297	z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
[1756]	298	END DO
[9019]	299	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
	300	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
[4990]	301	ENDIF
	302
	303	IF( iom_use("u_heattr") ) THEN
[9019]	304	z2d(:,:) = 0._wp
[4990]	305	DO jk = 1, jpkm1
	306	DO jj = 2, jpjm1
	307	DO ji = fs_2, fs_jpim1 ! vector opt.
	308	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
	309	END DO
	310	END DO
	311	END DO
[10425]	312	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9019]	313	CALL iom_put( "u_heattr", 0.5rcp z2d ) ! heat transport in i-direction
[4990]	314	ENDIF
[4761]	315
[4990]	316	IF( iom_use("u_salttr") ) THEN
[7753]	317	z2d(:,:) = 0.e0
[1756]	318	DO jk = 1, jpkm1
	319	DO jj = 2, jpjm1
	320	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	321	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
[1756]	322	END DO
	323	END DO
	324	END DO
[10425]	325	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9019]	326	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
[4990]	327	ENDIF
[4761]	328
[4990]	329
	330	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
[7753]	331	z3d(:,:,jpk) = 0.e0
[1756]	332	DO jk = 1, jpkm1
[7753]	333	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
[1756]	334	END DO
[9019]	335	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
[4990]	336	ENDIF
	337
	338	IF( iom_use("v_heattr") ) THEN
[7753]	339	z2d(:,:) = 0.e0
[4990]	340	DO jk = 1, jpkm1
	341	DO jj = 2, jpjm1
	342	DO ji = fs_2, fs_jpim1 ! vector opt.
	343	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
	344	END DO
	345	END DO
	346	END DO
[10425]	347	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9019]	348	CALL iom_put( "v_heattr", 0.5rcp z2d ) ! heat transport in j-direction
[4990]	349	ENDIF
[4761]	350
[4990]	351	IF( iom_use("v_salttr") ) THEN
[9019]	352	z2d(:,:) = 0._wp
[1756]	353	DO jk = 1, jpkm1
	354	DO jj = 2, jpjm1
	355	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	356	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
[1756]	357	END DO
	358	END DO
	359	END DO
[10425]	360	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9019]	361	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
[1756]	362	ENDIF
[7646]	363
	364	IF( iom_use("tosmint") ) THEN
[9019]	365	z2d(:,:) = 0._wp
[7646]	366	DO jk = 1, jpkm1
	367	DO jj = 2, jpjm1
	368	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	369	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
[7646]	370	END DO
	371	END DO
	372	END DO
[10425]	373	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9019]	374	CALL iom_put( "tosmint", rau0 * z2d ) ! Vertical integral of temperature
[7646]	375	ENDIF
	376	IF( iom_use("somint") ) THEN
[7753]	377	z2d(:,:)=0._wp
[7646]	378	DO jk = 1, jpkm1
	379	DO jj = 2, jpjm1
	380	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	381	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
[7646]	382	END DO
	383	END DO
	384	END DO
[10425]	385	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9019]	386	CALL iom_put( "somint", rau0 * z2d ) ! Vertical integral of salinity
[7646]	387	ENDIF
	388
[9019]	389	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
[2528]	390	!
[6140]	391
[9019]	392	IF (ln_diatmb) CALL dia_tmb ! tmb values
	393
	394	IF (ln_dia25h) CALL dia_25h( kt ) ! 25h averaging
[6140]	395
[9124]	396	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	397	!
[1482]	398	END SUBROUTINE dia_wri
	399
	400	#else
[2528]	401	!!----------------------------------------------------------------------
	402	!! Default option use IOIPSL library
	403	!!----------------------------------------------------------------------
	404
[9652]	405	INTEGER FUNCTION dia_wri_alloc()
	406	!!----------------------------------------------------------------------
	407	INTEGER, DIMENSION(2) :: ierr
	408	!!----------------------------------------------------------------------
	409	ierr = 0
	410	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	411	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	412	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
	413	!
	414	dia_wri_alloc = MAXVAL(ierr)
[10425]	415	CALL mpp_sum( 'diawri', dia_wri_alloc )
[9652]	416	!
	417	END FUNCTION dia_wri_alloc
	418
	419
[1561]	420	SUBROUTINE dia_wri( kt )
[3]	421	!!---------------------------------------------------------------------
	422	!! * ROUTINE dia_wri *
	423	!!
	424	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	425	!! NETCDF format is used by default
	426	!!
	427	!! ** Method : At the beginning of the first time step (nit000),
	428	!! define all the NETCDF files and fields
	429	!! At each time step call histdef to compute the mean if ncessary
	430	!! Each nwrite time step, output the instantaneous or mean fields
	431	!!----------------------------------------------------------------------
[5836]	432	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	433	!
[2528]	434	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	435	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	436	INTEGER :: inum = 11 ! temporary logical unit
[3294]	437	INTEGER :: ji, jj, jk ! dummy loop indices
	438	INTEGER :: ierr ! error code return from allocation
[2528]	439	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[3609]	440	INTEGER :: jn, ierror ! local integers
[6140]	441	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	442	!
[9019]	443	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
	444	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
[3]	445	!!----------------------------------------------------------------------
[3294]	446	!
[9124]	447	IF( ln_timing ) CALL timing_start('dia_wri')
[1482]	448	!
[9019]	449	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[10425]	450	CALL dia_wri_state( 'output.init' )
[1482]	451	ninist = 0
	452	ENDIF
	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
[6140]	463	zsto = nwrite * rdt
[1312]	464	clop = "inst("//TRIM(clop)//")"
[3]	465	#else
[6140]	466	zsto=rdt
[1312]	467	clop = "ave("//TRIM(clop)//")"
[3]	468	#endif
[6140]	469	zout = nwrite * rdt
	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
[895]	500	CALL dia_nam( clhstnam, nwrite,' ' )
[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
[3]	508	CALL dia_nam( clhstnam, nwrite, 'grid_T' )
	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
	544	CALL dia_nam( clhstnam, nwrite, 'grid_U' )
	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
	557	CALL dia_nam( clhstnam, nwrite, 'grid_V' ) ! filename
	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
	570	CALL dia_nam( clhstnam, nwrite, 'grid_W' ) ! filename
	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
[5407]	661	IF( .NOT. ln_cpl ) 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
[3]	669
[5407]	670	IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]	671	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	672	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	673	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	674	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	675	CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping" , "Kg/m2/s", & ! erp * sn
	676	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	677	ENDIF
	678
[3]	679	clmx ="l_max(only(x))" ! max index on a period
[5836]	680	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	681	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
[3]	682	#if defined key_diahth
	683	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	684	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	685	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	686	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	687	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	688	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[7646]	689	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
[3]	690	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	691	#endif
	692
[2528]	693	CALL histend( nid_T, snc4chunks=snc4set )
[3]	694
	695	! !!! nid_U : 3D
	696	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	697	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]	698	IF( ln_wave .AND. ln_sdw) THEN
	699	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	700	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	701	ENDIF
[3]	702	! !!! nid_U : 2D
[888]	703	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	704	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	705
[2528]	706	CALL histend( nid_U, snc4chunks=snc4set )
[3]	707
	708	! !!! nid_V : 3D
	709	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	710	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]	711	IF( ln_wave .AND. ln_sdw) THEN
	712	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	713	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	714	ENDIF
[3]	715	! !!! nid_V : 2D
[888]	716	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	717	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	718
[2528]	719	CALL histend( nid_V, snc4chunks=snc4set )
[3]	720
	721	! !!! nid_W : 3D
	722	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	723	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	724	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	725	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]	726	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
[255]	727	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	728
[9019]	729	IF( ln_zdfddm ) THEN
[3]	730	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	731	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	732	ENDIF
[7646]	733
	734	IF( ln_wave .AND. ln_sdw) THEN
	735	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	736	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	737	ENDIF
[3]	738	! !!! nid_W : 2D
[2528]	739	CALL histend( nid_W, snc4chunks=snc4set )
[3]	740
	741	IF(lwp) WRITE(numout,*)
	742	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	743	IF(ll_print) CALL FLUSH(numout )
	744
	745	ENDIF
	746
	747	! 2. Start writing data
	748	! ---------------------
	749
[4292]	750	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	751	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	752	! donne le nombre d'elements, et ndex la liste des indices a sortir
	753
[1334]	754	IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
[3]	755	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	756	WRITE(numout,*) '~~~~~~ '
	757	ENDIF
	758
[6140]	759	IF( .NOT.ln_linssh ) THEN
	760	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	761	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	762	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	763	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	764	ELSE
	765	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	766	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	767	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	768	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	769	ENDIF
[6140]	770	IF( .NOT.ln_linssh ) THEN
[7753]	771	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
[6140]	772	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	773	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
[4292]	774	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	775	ENDIF
[3]	776	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
[2528]	777	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	778	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	779	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	780	! (includes virtual salt flux beneath ice
	781	! in linear free surface case)
[6140]	782	IF( ln_linssh ) THEN
[7753]	783	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
[4292]	784	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
[7753]	785	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
[4292]	786	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	787	ENDIF
[888]	788	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	789	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	790	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	791	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	792	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	793	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[3609]	794	!
	795	IF( ln_icebergs ) THEN
	796	!
	797	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	798	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	799	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	800	!
	801	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	802	!
	803	IF( ln_bergdia ) THEN
	804	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	805	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	806	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	807	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	808	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	809	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	810	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	811	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	812	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	813	!
	814	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	815	ENDIF
	816	ENDIF
	817
[5407]	818	IF( .NOT. ln_cpl ) THEN
[4990]	819	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	820	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[7753]	821	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	822	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	823	ENDIF
[5407]	824	IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]	825	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	826	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[7753]	827	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	828	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	829	ENDIF
	830	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	831	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	832
	833	#if defined key_diahth
	834	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	835	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	836	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	837	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	838	#endif
[888]	839
[3]	840	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
[888]	841	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	842
	843	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
[888]	844	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	845
	846	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
	847	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	848	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	849	IF( ln_zdfddm ) THEN
	850	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	851	ENDIF
	852
[7646]	853	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	854	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	855	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	856	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	857	ENDIF
	858
[1318]	859	! 3. Close all files
	860	! ---------------------------------------
[1561]	861	IF( kt == nitend ) THEN
[3]	862	CALL histclo( nid_T )
	863	CALL histclo( nid_U )
	864	CALL histclo( nid_V )
	865	CALL histclo( nid_W )
	866	ENDIF
[2528]	867	!
[9124]	868	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	869	!
[3]	870	END SUBROUTINE dia_wri
[1567]	871	#endif
	872
[10425]	873	SUBROUTINE dia_wri_state( cdfile_name )
[3]	874	!!---------------------------------------------------------------------
	875	!! * ROUTINE dia_wri_state *
	876	!!
	877	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	878	!! the instantaneous ocean state and forcing fields.
	879	!! Used to find errors in the initial state or save the last
	880	!! ocean state in case of abnormal end of a simulation
	881	!!
	882	!! ** Method : NetCDF files using ioipsl
	883	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	884	!! File 'output.abort.nc' is created in case of abnormal job end
	885	!!----------------------------------------------------------------------
[1334]	886	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[11521]	887	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zisfdebug
[10425]	888	!!
[11521]	889	INTEGER :: inum, jk
[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
[11876]	909	CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep ) ! now k-velocity
	910	CALL iom_rstput( 0, 0, inum, 'ht_n' , ht_n ) ! now k-velocity
[11423]	911	IF ( ln_isf ) THEN
	912	IF (ln_isfcav_mlt) THEN
	913	CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav ) ! now k-velocity
	914	CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav ) ! now k-velocity
	915	CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav ) ! now k-velocity
	916	CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8) ) ! now k-velocity
	917	CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8) ) ! now k-velocity
	918	END IF
	919	IF (ln_isfpar_mlt) THEN
[11494]	920	CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,8) ) ! now k-velocity
[11423]	921	CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par ) ! now k-velocity
	922	CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par ) ! now k-velocity
	923	CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par ) ! now k-velocity
	924	CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8) ) ! now k-velocity
	925	CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8) ) ! now k-velocity
	926	END IF
	927	END IF
[11876]	928
	929	IF ( ln_isf ) THEN
	930	IF (ln_isfcav_mlt) CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 )
	931	IF (ln_isfpar_mlt) CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 )
	932	END IF
	933
[7646]	934	IF( ALLOCATED(ahtu) ) THEN
[10425]	935	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	936	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	937	ENDIF
	938	IF( ALLOCATED(ahmt) ) THEN
[10425]	939	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	940	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	941	ENDIF
[10425]	942	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	943	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	944	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	945	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	946	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	947	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[6140]	948	IF( .NOT.ln_linssh ) THEN
[10425]	949	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
	950	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
	951	END IF
[7646]	952	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	953	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	954	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	955	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	956	ENDIF
[10425]	957
	958	#if defined key_si3
	959	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
	960	CALL ice_wri_state( inum )
[1561]	961	ENDIF
	962	#endif
[10425]	963	!
	964	CALL iom_close( inum )
[3294]	965	!
[3]	966	END SUBROUTINE dia_wri_state
[6140]	967
[3]	968	!!======================================================================
	969	END MODULE diawri

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

Download in other formats: