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 NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DIA/diawri.F90 @ 13085

Last change on this file since 13085 was 13085, checked in by davestorkey, 4 years ago
r4.0-HEAD : fix for ticket #2425.
Property svn:keywords set to `Id`
File size: 49.2 KB

Rev	Line
[3]	1	MODULE diawri
	2	!!======================================================================
	3	!! * MODULE diawri *
	4	!! Ocean diagnostics : write ocean output files
	5	!!=====================================================================
[2528]	6	!! History : OPA ! 1991-03 (M.-A. Foujols) Original code
	7	!! 4.0 ! 1991-11 (G. Madec)
	8	!! ! 1992-06 (M. Imbard) correction restart file
	9	!! ! 1992-07 (M. Imbard) split into diawri and rstwri
	10	!! ! 1993-03 (M. Imbard) suppress writibm
	11	!! ! 1998-01 (C. Levy) NETCDF format using ioipsl INTERFACE
	12	!! ! 1999-02 (E. Guilyardi) name of netCDF files + variables
	13	!! 8.2 ! 2000-06 (M. Imbard) Original code (diabort.F)
	14	!! NEMO 1.0 ! 2002-06 (A.Bozec, E. Durand) Original code (diainit.F)
	15	!! - ! 2002-09 (G. Madec) F90: Free form and module
	16	!! - ! 2002-12 (G. Madec) merge of diabort and diainit, F90
	17	!! ! 2005-11 (V. Garnier) Surface pressure gradient organization
	18	!! 3.2 ! 2008-11 (B. Lemaire) creation from old diawri
[5836]	19	!! 3.7 ! 2014-01 (G. Madec) remove eddy induced velocity from no-IOM output
	20	!! ! change name of output variables in dia_wri_state
[2528]	21	!!----------------------------------------------------------------------
[3]	22
	23	!!----------------------------------------------------------------------
[2528]	24	!! dia_wri : create the standart output files
	25	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
	26	!!----------------------------------------------------------------------
[9019]	27	USE oce ! ocean dynamics and tracers
	28	USE dom_oce ! ocean space and time domain
	29	USE phycst ! physical constants
	30	USE dianam ! build name of file (routine)
	31	USE diahth ! thermocline diagnostics
	32	USE dynadv , ONLY: ln_dynadv_vec
	33	USE icb_oce ! Icebergs
	34	USE icbdia ! Iceberg budgets
	35	USE ldftra ! lateral physics: eddy diffusivity coef.
	36	USE ldfdyn ! lateral physics: eddy viscosity coef.
	37	USE sbc_oce ! Surface boundary condition: ocean fields
	38	USE sbc_ice ! Surface boundary condition: ice fields
	39	USE sbcssr ! restoring term toward SST/SSS climatology
	40	USE sbcwave ! wave parameters
	41	USE wet_dry ! wetting and drying
	42	USE zdf_oce ! ocean vertical physics
	43	USE zdfdrg ! ocean vertical physics: top/bottom friction
	44	USE zdfmxl ! mixed layer
[6140]	45	!
[9019]	46	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	47	USE in_out_manager ! I/O manager
	48	USE dia25h ! 25h Mean output
	49	USE iom !
	50	USE ioipsl !
[5463]	51
[9570]	52	#if defined key_si3
[10425]	53	USE ice
[9019]	54	USE icewri
[1482]	55	#endif
[2715]	56	USE lib_mpp ! MPP library
[3294]	57	USE timing ! preformance summary
[6140]	58	USE diurnal_bulk ! diurnal warm layer
	59	USE cool_skin ! Cool skin
[2528]	60
[3]	61	IMPLICIT NONE
	62	PRIVATE
	63
[2528]	64	PUBLIC dia_wri ! routines called by step.F90
	65	PUBLIC dia_wri_state
[2715]	66	PUBLIC dia_wri_alloc ! Called by nemogcm module
[3]	67
[2528]	68	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[3609]	69	INTEGER :: nb_T , ndim_bT ! grid_T file
[2528]	70	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	71	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	72	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
	73	INTEGER :: ndex(1) ! ???
[2715]	74	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
	75	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]	76	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	77
	78	!! * Substitutions
[1756]	79	# include "vectopt_loop_substitute.h90"
[3]	80	!!----------------------------------------------------------------------
[9598]	81	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	82	!! $Id$
[10068]	83	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	84	!!----------------------------------------------------------------------
	85	CONTAINS
	86
[6140]	87	#if defined key_iomput
[3]	88	!!----------------------------------------------------------------------
[2528]	89	!! 'key_iomput' use IOM library
	90	!!----------------------------------------------------------------------
[9652]	91	INTEGER FUNCTION dia_wri_alloc()
	92	!
	93	dia_wri_alloc = 0
	94	!
	95	END FUNCTION dia_wri_alloc
[2715]	96
[9652]	97
[1561]	98	SUBROUTINE dia_wri( kt )
[1482]	99	!!---------------------------------------------------------------------
	100	!! * ROUTINE dia_wri *
	101	!!
	102	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	103	!! NETCDF format is used by default
	104	!!
	105	!! ** Method : use iom_put
	106	!!----------------------------------------------------------------------
	107	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1756]	108	!!
[9019]	109	INTEGER :: ji, jj, jk ! dummy loop indices
	110	INTEGER :: ikbot ! local integer
	111	REAL(wp):: zztmp , zztmpx ! local scalar
	112	REAL(wp):: zztmp2, zztmpy ! - -
	113	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
	114	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
[1482]	115	!!----------------------------------------------------------------------
	116	!
[9124]	117	IF( ln_timing ) CALL timing_start('dia_wri')
[3294]	118	!
[1482]	119	! Output the initial state and forcings
	120	IF( ninist == 1 ) THEN
[10425]	121	CALL dia_wri_state( 'output.init' )
[1482]	122	ninist = 0
	123	ENDIF
[3]	124
[6351]	125	! Output of initial vertical scale factor
	126	CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10114]	127	CALL iom_put("e3u_0", e3u_0(:,:,:) )
	128	CALL iom_put("e3v_0", e3v_0(:,:,:) )
[6351]	129	!
[6387]	130	CALL iom_put( "e3t" , e3t_n(:,:,:) )
	131	CALL iom_put( "e3u" , e3u_n(:,:,:) )
	132	CALL iom_put( "e3v" , e3v_n(:,:,:) )
	133	CALL iom_put( "e3w" , e3w_n(:,:,:) )
[6351]	134	IF( iom_use("e3tdef") ) &
[6387]	135	CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
[5461]	136
[9023]	137	IF( ll_wd ) THEN
	138	CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)
	139	ELSE
	140	CALL iom_put( "ssh" , sshn ) ! sea surface height
	141	ENDIF
	142
[7646]	143	IF( iom_use("wetdep") ) & ! wet depth
[9023]	144	CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
[5107]	145
	146	CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature
	147	CALL iom_put( "sst", tsn(:,:,1,jp_tem) ) ! surface temperature
	148	IF ( iom_use("sbt") ) THEN
[4990]	149	DO jj = 1, jpj
	150	DO ji = 1, jpi
[9019]	151	ikbot = mbkt(ji,jj)
	152	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
[4990]	153	END DO
[5107]	154	END DO
	155	CALL iom_put( "sbt", z2d ) ! bottom temperature
	156	ENDIF
	157
	158	CALL iom_put( "soce", tsn(:,:,:,jp_sal) ) ! 3D salinity
	159	CALL iom_put( "sss", tsn(:,:,1,jp_sal) ) ! surface salinity
	160	IF ( iom_use("sbs") ) THEN
[4990]	161	DO jj = 1, jpj
	162	DO ji = 1, jpi
[9019]	163	ikbot = mbkt(ji,jj)
	164	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
[4990]	165	END DO
[5107]	166	END DO
	167	CALL iom_put( "sbs", z2d ) ! bottom salinity
	168	ENDIF
[5463]	169
[13085]	170	CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
	171
[5463]	172	IF ( iom_use("taubot") ) THEN ! bottom stress
[9019]	173	zztmp = rau0 * 0.25
[7753]	174	z2d(:,:) = 0._wp
[5463]	175	DO jj = 2, jpjm1
	176	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	177	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * un(ji ,jj,mbku(ji ,jj)) )**2 &
	178	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj)) )**2 &
	179	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vn(ji,jj ,mbkv(ji,jj )) )**2 &
	180	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1)) )**2
	181	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
[5463]	182	!
[9019]	183	END DO
	184	END DO
[10425]	185	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[5463]	186	CALL iom_put( "taubot", z2d )
	187	ENDIF
[5107]	188
[9019]	189	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
	190	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
[5107]	191	IF ( iom_use("sbu") ) THEN
[4990]	192	DO jj = 1, jpj
	193	DO ji = 1, jpi
[9019]	194	ikbot = mbku(ji,jj)
	195	z2d(ji,jj) = un(ji,jj,ikbot)
[4990]	196	END DO
[5107]	197	END DO
	198	CALL iom_put( "sbu", z2d ) ! bottom i-current
	199	ENDIF
	200
[9019]	201	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
	202	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
[5107]	203	IF ( iom_use("sbv") ) THEN
[4990]	204	DO jj = 1, jpj
	205	DO ji = 1, jpi
[9019]	206	ikbot = mbkv(ji,jj)
	207	z2d(ji,jj) = vn(ji,jj,ikbot)
[4990]	208	END DO
[5107]	209	END DO
	210	CALL iom_put( "sbv", z2d ) ! bottom j-current
[4990]	211	ENDIF
[1482]	212
[11418]	213	IF( ln_zad_Aimp ) wn = wn + wi ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
	214	!
[5461]	215	CALL iom_put( "woce", wn ) ! vertical velocity
	216	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
	217	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[7753]	218	z2d(:,:) = rau0 * e1e2t(:,:)
[5461]	219	DO jk = 1, jpk
[7753]	220	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
[5461]	221	END DO
	222	CALL iom_put( "w_masstr" , z3d )
	223	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
	224	ENDIF
[11418]	225	!
	226	IF( ln_zad_Aimp ) wn = wn - wi ! Remove implicit part of vertical velocity that was added for diagnostic output
[5461]	227
[9019]	228	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	229	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
	230	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
[5107]	231
[9019]	232	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
	233	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
[6351]	234
[5107]	235	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
[4990]	236	DO jj = 2, jpjm1 ! sst gradient
	237	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	238	zztmp = tsn(ji,jj,1,jp_tem)
	239	zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,1,jp_tem) ) * r1_e1u(ji-1,jj)
	240	zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
[4990]	241	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	242	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
	243	END DO
[1756]	244	END DO
[10425]	245	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[9019]	246	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
[7753]	247	z2d(:,:) = SQRT( z2d(:,:) )
[9019]	248	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
[4990]	249	ENDIF
	250
[9019]	251	! heat and salt contents
[4990]	252	IF( iom_use("heatc") ) THEN
[7753]	253	z2d(:,:) = 0._wp
[4990]	254	DO jk = 1, jpkm1
[5107]	255	DO jj = 1, jpj
	256	DO ji = 1, jpi
[6140]	257	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
[4990]	258	END DO
[4761]	259	END DO
	260	END DO
[9019]	261	CALL iom_put( "heatc", rau0_rcp * z2d ) ! vertically integrated heat content (J/m2)
[4990]	262	ENDIF
	263
	264	IF( iom_use("saltc") ) THEN
[7753]	265	z2d(:,:) = 0._wp
[4990]	266	DO jk = 1, jpkm1
[5107]	267	DO jj = 1, jpj
	268	DO ji = 1, jpi
[6140]	269	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
[4990]	270	END DO
	271	END DO
	272	END DO
[9019]	273	CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
[4990]	274	ENDIF
[4840]	275	!
[4990]	276	IF ( iom_use("eken") ) THEN
[9399]	277	z3d(:,:,jpk) = 0._wp
[4990]	278	DO jk = 1, jpkm1
	279	DO jj = 2, jpjm1
	280	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	281	zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
	282	z3d(ji,jj,jk) = zztmp * ( un(ji-1,jj,jk)*2 e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) &
	283	& + un(ji ,jj,jk)*2 e2u(ji ,jj) * e3u_n(ji ,jj,jk) &
	284	& + vn(ji,jj-1,jk)*2 e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) &
	285	& + vn(ji,jj ,jk)*2 e1v(ji,jj ) * e3v_n(ji,jj ,jk) )
	286	END DO
	287	END DO
	288	END DO
[10425]	289	CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
[9019]	290	CALL iom_put( "eken", z3d ) ! kinetic energy
[4990]	291	ENDIF
[6351]	292	!
	293	CALL iom_put( "hdiv", hdivn ) ! Horizontal divergence
	294	!
[7646]	295	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
[7753]	296	z3d(:,:,jpk) = 0.e0
	297	z2d(:,:) = 0.e0
[1756]	298	DO jk = 1, jpkm1
[7753]	299	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
	300	z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
[1756]	301	END DO
[9019]	302	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
	303	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
[4990]	304	ENDIF
	305
	306	IF( iom_use("u_heattr") ) THEN
[9019]	307	z2d(:,:) = 0._wp
[4990]	308	DO jk = 1, jpkm1
	309	DO jj = 2, jpjm1
	310	DO ji = fs_2, fs_jpim1 ! vector opt.
	311	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
	312	END DO
	313	END DO
	314	END DO
[10425]	315	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9019]	316	CALL iom_put( "u_heattr", 0.5rcp z2d ) ! heat transport in i-direction
[4990]	317	ENDIF
[4761]	318
[4990]	319	IF( iom_use("u_salttr") ) THEN
[7753]	320	z2d(:,:) = 0.e0
[1756]	321	DO jk = 1, jpkm1
	322	DO jj = 2, jpjm1
	323	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	324	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
[1756]	325	END DO
	326	END DO
	327	END DO
[10425]	328	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
[9019]	329	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
[4990]	330	ENDIF
[4761]	331
[4990]	332
	333	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
[7753]	334	z3d(:,:,jpk) = 0.e0
[1756]	335	DO jk = 1, jpkm1
[7753]	336	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
[1756]	337	END DO
[9019]	338	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
[4990]	339	ENDIF
	340
	341	IF( iom_use("v_heattr") ) THEN
[7753]	342	z2d(:,:) = 0.e0
[4990]	343	DO jk = 1, jpkm1
	344	DO jj = 2, jpjm1
	345	DO ji = fs_2, fs_jpim1 ! vector opt.
	346	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
	347	END DO
	348	END DO
	349	END DO
[10425]	350	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9019]	351	CALL iom_put( "v_heattr", 0.5rcp z2d ) ! heat transport in j-direction
[4990]	352	ENDIF
[4761]	353
[4990]	354	IF( iom_use("v_salttr") ) THEN
[9019]	355	z2d(:,:) = 0._wp
[1756]	356	DO jk = 1, jpkm1
	357	DO jj = 2, jpjm1
	358	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	359	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
[1756]	360	END DO
	361	END DO
	362	END DO
[10425]	363	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
[9019]	364	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
[1756]	365	ENDIF
[7646]	366
	367	IF( iom_use("tosmint") ) THEN
[9019]	368	z2d(:,:) = 0._wp
[7646]	369	DO jk = 1, jpkm1
	370	DO jj = 2, jpjm1
	371	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	372	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
[7646]	373	END DO
	374	END DO
	375	END DO
[10425]	376	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9019]	377	CALL iom_put( "tosmint", rau0 * z2d ) ! Vertical integral of temperature
[7646]	378	ENDIF
	379	IF( iom_use("somint") ) THEN
[7753]	380	z2d(:,:)=0._wp
[7646]	381	DO jk = 1, jpkm1
	382	DO jj = 2, jpjm1
	383	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	384	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
[7646]	385	END DO
	386	END DO
	387	END DO
[10425]	388	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
[9019]	389	CALL iom_put( "somint", rau0 * z2d ) ! Vertical integral of salinity
[7646]	390	ENDIF
	391
[9019]	392	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
[2528]	393	!
[9019]	394
	395	IF (ln_dia25h) CALL dia_25h( kt ) ! 25h averaging
[6140]	396
[9124]	397	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	398	!
[1482]	399	END SUBROUTINE dia_wri
	400
	401	#else
[2528]	402	!!----------------------------------------------------------------------
	403	!! Default option use IOIPSL library
	404	!!----------------------------------------------------------------------
	405
[9652]	406	INTEGER FUNCTION dia_wri_alloc()
	407	!!----------------------------------------------------------------------
	408	INTEGER, DIMENSION(2) :: ierr
	409	!!----------------------------------------------------------------------
[12494]	410	IF( nn_write == -1 ) THEN
	411	dia_wri_alloc = 0
	412	ELSE
	413	ierr = 0
	414	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	415	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	416	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
[9652]	417	!
[12494]	418	dia_wri_alloc = MAXVAL(ierr)
	419	CALL mpp_sum( 'diawri', dia_wri_alloc )
	420	!
	421	ENDIF
[9652]	422	!
	423	END FUNCTION dia_wri_alloc
	424
	425
[1561]	426	SUBROUTINE dia_wri( kt )
[3]	427	!!---------------------------------------------------------------------
	428	!! * ROUTINE dia_wri *
	429	!!
	430	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	431	!! NETCDF format is used by default
	432	!!
	433	!! ** Method : At the beginning of the first time step (nit000),
	434	!! define all the NETCDF files and fields
	435	!! At each time step call histdef to compute the mean if ncessary
[11536]	436	!! Each nn_write time step, output the instantaneous or mean fields
[3]	437	!!----------------------------------------------------------------------
[5836]	438	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	439	!
[2528]	440	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	441	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	442	INTEGER :: inum = 11 ! temporary logical unit
[3294]	443	INTEGER :: ji, jj, jk ! dummy loop indices
	444	INTEGER :: ierr ! error code return from allocation
[2528]	445	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[3609]	446	INTEGER :: jn, ierror ! local integers
[6140]	447	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	448	!
[9019]	449	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
	450	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
[3]	451	!!----------------------------------------------------------------------
[1482]	452	!
[9019]	453	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[10425]	454	CALL dia_wri_state( 'output.init' )
[1482]	455	ninist = 0
	456	ENDIF
	457	!
[11536]	458	IF( nn_write == -1 ) RETURN ! we will never do any output
	459	!
	460	IF( ln_timing ) CALL timing_start('dia_wri')
	461	!
[3]	462	! 0. Initialisation
	463	! -----------------
[632]	464
[9019]	465	ll_print = .FALSE. ! local variable for debugging
[3]	466	ll_print = ll_print .AND. lwp
	467
	468	! Define frequency of output and means
[5566]	469	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]	470	#if defined key_diainstant
[11536]	471	zsto = nn_write * rdt
[1312]	472	clop = "inst("//TRIM(clop)//")"
[3]	473	#else
[6140]	474	zsto=rdt
[1312]	475	clop = "ave("//TRIM(clop)//")"
[3]	476	#endif
[11536]	477	zout = nn_write * rdt
[6140]	478	zmax = ( nitend - nit000 + 1 ) * rdt
[3]	479
	480	! Define indices of the horizontal output zoom and vertical limit storage
	481	iimi = 1 ; iima = jpi
	482	ijmi = 1 ; ijma = jpj
	483	ipk = jpk
	484
	485	! define time axis
[1334]	486	it = kt
	487	itmod = kt - nit000 + 1
[3]	488
	489
	490	! 1. Define NETCDF files and fields at beginning of first time step
	491	! -----------------------------------------------------------------
	492
	493	IF( kt == nit000 ) THEN
	494
	495	! Define the NETCDF files (one per grid)
[632]	496
[3]	497	! Compute julian date from starting date of the run
[1309]	498	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	499	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	500	IF(lwp)WRITE(numout,*)
	501	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	502	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	503	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	504	' limit storage in depth = ', ipk
	505
	506	! WRITE root name in date.file for use by postpro
[1581]	507	IF(lwp) THEN
[11536]	508	CALL dia_nam( clhstnam, nn_write,' ' )
[1581]	509	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	510	WRITE(inum,*) clhstnam
	511	CLOSE(inum)
	512	ENDIF
[632]	513
[3]	514	! Define the T grid FILE ( nid_T )
[632]	515
[11536]	516	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[3]	517	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	518	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	519	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	520	& nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	521	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	522	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	523	! ! Index of ocean points
	524	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	525	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	526	!
	527	IF( ln_icebergs ) THEN
	528	!
	529	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	530	!! that routine is called from nemogcm, so do it here immediately before its needed
	531	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	532	CALL mpp_sum( 'diawri', ierror )
[3609]	533	IF( ierror /= 0 ) THEN
	534	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	535	RETURN
	536	ENDIF
	537	!
	538	!! iceberg vertical coordinate is class number
	539	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	540	& "number", nclasses, class_num, nb_T )
	541	!
	542	!! each class just needs the surface index pattern
	543	ndim_bT = 3
	544	DO jn = 1,nclasses
	545	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	546	ENDDO
	547	!
	548	ENDIF
[3]	549
	550	! Define the U grid FILE ( nid_U )
	551
[11536]	552	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[3]	553	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	554	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	555	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	556	& nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	557	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	558	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	559	! ! Index of ocean points
	560	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	561	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	562
	563	! Define the V grid FILE ( nid_V )
	564
[11536]	565	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[3]	566	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	567	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	568	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	569	& nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	570	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	571	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	572	! ! Index of ocean points
	573	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	574	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	575
	576	! Define the W grid FILE ( nid_W )
	577
[11536]	578	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[3]	579	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	580	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	581	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	582	& nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	583	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	584	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	585
[632]	586
[3]	587	! Declare all the output fields as NETCDF variables
	588
	589	! !!! nid_T : 3D
	590	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	591	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	592	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	593	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]	594	IF( .NOT.ln_linssh ) THEN
[4292]	595	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	596	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	597	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	598	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	599	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	600	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	601	ENDIF
[3]	602	! !!! nid_T : 2D
	603	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	604	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	605	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	606	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	607	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	608	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	609	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	610	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	611	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	612	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	613	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	614	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[6140]	615	IF( ln_linssh ) THEN
[4292]	616	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
[3625]	617	& , "KgC/m2/s", & ! sosst_cd
[4292]	618	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	619	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
[3625]	620	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	621	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	622	ENDIF
[888]	623	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	624	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	625	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	626	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	627	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	628	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	629	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	630	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	631	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	632	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	633	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	634	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3609]	635	!
	636	IF( ln_icebergs ) THEN
	637	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	638	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	639	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	640	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	641	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	642	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	643	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	644	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	645	IF( ln_bergdia ) THEN
	646	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	647	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	648	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	649	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	650	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	651	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	652	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	653	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	654	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	655	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	656	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	657	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	658	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	659	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	660	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	661	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	662	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	663	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	664	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	665	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	666	ENDIF
	667	ENDIF
	668
[11536]	669	IF( ln_ssr ) THEN
[4990]	670	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	671	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	672	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	673	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	674	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	675	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	676	ENDIF
[11536]	677
[3]	678	clmx ="l_max(only(x))" ! max index on a period
[5836]	679	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	680	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
[3]	681	#if defined key_diahth
	682	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	683	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	684	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	685	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	686	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	687	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[7646]	688	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
[3]	689	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	690	#endif
	691
[2528]	692	CALL histend( nid_T, snc4chunks=snc4set )
[3]	693
	694	! !!! nid_U : 3D
	695	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	696	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]	697	IF( ln_wave .AND. ln_sdw) THEN
	698	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	699	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	700	ENDIF
[3]	701	! !!! nid_U : 2D
[888]	702	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	703	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	704
[2528]	705	CALL histend( nid_U, snc4chunks=snc4set )
[3]	706
	707	! !!! nid_V : 3D
	708	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	709	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]	710	IF( ln_wave .AND. ln_sdw) THEN
	711	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	712	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	713	ENDIF
[3]	714	! !!! nid_V : 2D
[888]	715	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	716	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	717
[2528]	718	CALL histend( nid_V, snc4chunks=snc4set )
[3]	719
	720	! !!! nid_W : 3D
	721	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	722	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	723	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	724	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]	725	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
[255]	726	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	727
[9019]	728	IF( ln_zdfddm ) THEN
[3]	729	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	730	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	731	ENDIF
[7646]	732
	733	IF( ln_wave .AND. ln_sdw) THEN
	734	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	735	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	736	ENDIF
[3]	737	! !!! nid_W : 2D
[2528]	738	CALL histend( nid_W, snc4chunks=snc4set )
[3]	739
	740	IF(lwp) WRITE(numout,*)
	741	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	742	IF(ll_print) CALL FLUSH(numout )
	743
	744	ENDIF
	745
	746	! 2. Start writing data
	747	! ---------------------
	748
[4292]	749	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	750	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	751	! donne le nombre d'elements, et ndex la liste des indices a sortir
	752
[11536]	753	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
[3]	754	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	755	WRITE(numout,*) '~~~~~~ '
	756	ENDIF
	757
[6140]	758	IF( .NOT.ln_linssh ) THEN
	759	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	760	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	761	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	762	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	763	ELSE
	764	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	765	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	766	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	767	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	768	ENDIF
[6140]	769	IF( .NOT.ln_linssh ) THEN
[7753]	770	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
[6140]	771	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	772	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
[4292]	773	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	774	ENDIF
[3]	775	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
[2528]	776	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	777	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	778	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	779	! (includes virtual salt flux beneath ice
	780	! in linear free surface case)
[6140]	781	IF( ln_linssh ) THEN
[7753]	782	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
[4292]	783	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
[7753]	784	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
[4292]	785	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	786	ENDIF
[888]	787	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	788	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	789	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	790	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	791	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	792	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[3609]	793	!
	794	IF( ln_icebergs ) THEN
	795	!
	796	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	797	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	798	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	799	!
	800	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	801	!
	802	IF( ln_bergdia ) THEN
	803	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	804	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	805	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	806	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	807	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	808	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	809	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	810	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	811	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	812	!
	813	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	814	ENDIF
	815	ENDIF
	816
[11536]	817	IF( ln_ssr ) THEN
[4990]	818	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	819	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[11536]	820	zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	821	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	822	ENDIF
	823	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	824	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	825
	826	#if defined key_diahth
	827	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	828	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	829	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	830	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	831	#endif
[888]	832
[3]	833	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
[888]	834	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	835
	836	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
[888]	837	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	838
[11418]	839	IF( ln_zad_Aimp ) THEN
	840	CALL histwrite( nid_W, "vovecrtz", it, wn + wi , ndim_T, ndex_T ) ! vert. current
	841	ELSE
	842	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
	843	ENDIF
[3]	844	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	845	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	846	IF( ln_zdfddm ) THEN
	847	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	848	ENDIF
	849
[7646]	850	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	851	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	852	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	853	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	854	ENDIF
	855
[1318]	856	! 3. Close all files
	857	! ---------------------------------------
[1561]	858	IF( kt == nitend ) THEN
[3]	859	CALL histclo( nid_T )
	860	CALL histclo( nid_U )
	861	CALL histclo( nid_V )
	862	CALL histclo( nid_W )
	863	ENDIF
[2528]	864	!
[9124]	865	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	866	!
[3]	867	END SUBROUTINE dia_wri
[1567]	868	#endif
	869
[10425]	870	SUBROUTINE dia_wri_state( cdfile_name )
[3]	871	!!---------------------------------------------------------------------
	872	!! * ROUTINE dia_wri_state *
	873	!!
	874	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	875	!! the instantaneous ocean state and forcing fields.
	876	!! Used to find errors in the initial state or save the last
	877	!! ocean state in case of abnormal end of a simulation
	878	!!
	879	!! ** Method : NetCDF files using ioipsl
	880	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	881	!! File 'output.abort.nc' is created in case of abnormal job end
	882	!!----------------------------------------------------------------------
[1334]	883	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	884	!!
	885	INTEGER :: inum
[3]	886	!!----------------------------------------------------------------------
[3294]	887	!
[648]	888	IF(lwp) WRITE(numout,*)
	889	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	890	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
[10425]	891	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
[648]	892
[9570]	893	#if defined key_si3
[10425]	894	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
[1482]	895	#else
[10425]	896	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
[1482]	897	#endif
[3]	898
[10425]	899	CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) ) ! now temperature
	900	CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) ) ! now salinity
	901	CALL iom_rstput( 0, 0, inum, 'sossheig', sshn ) ! sea surface height
	902	CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity
	903	CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity
[11418]	904	IF( ln_zad_Aimp ) THEN
	905	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi ) ! now k-velocity
	906	ELSE
	907	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity
	908	ENDIF
[7646]	909	IF( ALLOCATED(ahtu) ) THEN
[10425]	910	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	911	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	912	ENDIF
	913	IF( ALLOCATED(ahmt) ) THEN
[10425]	914	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	915	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	916	ENDIF
[10425]	917	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	918	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	919	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	920	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	921	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	922	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[6140]	923	IF( .NOT.ln_linssh ) THEN
[10425]	924	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
	925	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
	926	END IF
[7646]	927	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	928	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	929	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	930	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	931	ENDIF
[10425]	932
	933	#if defined key_si3
	934	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
	935	CALL ice_wri_state( inum )
[1561]	936	ENDIF
	937	#endif
[10425]	938	!
	939	CALL iom_close( inum )
[3294]	940	!
[3]	941	END SUBROUTINE dia_wri_state
[6140]	942
[3]	943	!!======================================================================
	944	END MODULE diawri

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

Download in other formats: