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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5217

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

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

Download in other formats: