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 @ 3868

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

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

Download in other formats: