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

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

diawri.F90 in branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5098

Last change on this file since 5098 was 5098, checked in by mathiot, 9 years ago
add wmask, wumask, wvmask and restore loop order and add flag to ignore specific isf code if no isf
Property svn:keywords set to `Id`
File size: 54.4 KB

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

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

Download in other formats: