Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

diawri.F90 @ 11533

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

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

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

Context Navigation

source: NEMO/releases/release-4.0/src/OCE/DIA/diawri.F90 @ 11533

Download in other formats: