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 branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5866

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

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

Download in other formats: