New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diawri.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4378

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

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

Download in other formats: