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_CMCC_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6051

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

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

Download in other formats: