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

Context Navigation

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6010

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

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

Download in other formats: