Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

diawri.F90 @ 12386

Last change on this file since 12386 was 12386, checked in by cguiavarch, 4 years ago
corrections to tramle & zdfosm to avoid division by 0 and restored OSMOSIS variables in output.abort
File size: 50.0 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
[12386]	45	USE zdfosm ! mixed layer
[6140]	46	!
[9019]	47	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	48	USE in_out_manager ! I/O manager
	49	USE diatmb ! Top,middle,bottom output
	50	USE dia25h ! 25h Mean output
	51	USE iom !
	52	USE ioipsl !
[5463]	53
[9570]	54	#if defined key_si3
[10425]	55	USE ice
[9019]	56	USE icewri
[1482]	57	#endif
[2715]	58	USE lib_mpp ! MPP library
[3294]	59	USE timing ! preformance summary
[6140]	60	USE diurnal_bulk ! diurnal warm layer
	61	USE cool_skin ! Cool skin
[2528]	62
[3]	63	IMPLICIT NONE
	64	PRIVATE
	65
[2528]	66	PUBLIC dia_wri ! routines called by step.F90
	67	PUBLIC dia_wri_state
[2715]	68	PUBLIC dia_wri_alloc ! Called by nemogcm module
[3]	69
[2528]	70	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[3609]	71	INTEGER :: nb_T , ndim_bT ! grid_T file
[2528]	72	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	73	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	74	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
	75	INTEGER :: ndex(1) ! ???
[2715]	76	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
	77	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]	78	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	79
	80	!! * Substitutions
[1756]	81	# include "vectopt_loop_substitute.h90"
[3]	82	!!----------------------------------------------------------------------
[9598]	83	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	84	!! $Id$
[10068]	85	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	86	!!----------------------------------------------------------------------
	87	CONTAINS
	88
[6140]	89	#if defined key_iomput
[3]	90	!!----------------------------------------------------------------------
[2528]	91	!! 'key_iomput' use IOM library
	92	!!----------------------------------------------------------------------
[9652]	93	INTEGER FUNCTION dia_wri_alloc()
	94	!
	95	dia_wri_alloc = 0
	96	!
	97	END FUNCTION dia_wri_alloc
[2715]	98
[9652]	99
[1561]	100	SUBROUTINE dia_wri( kt )
[1482]	101	!!---------------------------------------------------------------------
	102	!! * ROUTINE dia_wri *
	103	!!
	104	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	105	!! NETCDF format is used by default
	106	!!
	107	!! ** Method : use iom_put
	108	!!----------------------------------------------------------------------
	109	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1756]	110	!!
[9019]	111	INTEGER :: ji, jj, jk ! dummy loop indices
	112	INTEGER :: ikbot ! local integer
	113	REAL(wp):: zztmp , zztmpx ! local scalar
	114	REAL(wp):: zztmp2, zztmpy ! - -
	115	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
	116	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
[1482]	117	!!----------------------------------------------------------------------
	118	!
[9124]	119	IF( ln_timing ) CALL timing_start('dia_wri')
[3294]	120	!
[1482]	121	! Output the initial state and forcings
	122	IF( ninist == 1 ) THEN
[10425]	123	CALL dia_wri_state( 'output.init' )
[1482]	124	ninist = 0
	125	ENDIF
[3]	126
[6351]	127	! Output of initial vertical scale factor
	128	CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10114]	129	CALL iom_put("e3u_0", e3u_0(:,:,:) )
	130	CALL iom_put("e3v_0", e3v_0(:,:,:) )
[6351]	131	!
[6387]	132	CALL iom_put( "e3t" , e3t_n(:,:,:) )
	133	CALL iom_put( "e3u" , e3u_n(:,:,:) )
	134	CALL iom_put( "e3v" , e3v_n(:,:,:) )
	135	CALL iom_put( "e3w" , e3w_n(:,:,:) )
[6351]	136	IF( iom_use("e3tdef") ) &
[6387]	137	CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
[5461]	138
[9023]	139	IF( ll_wd ) THEN
	140	CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)
	141	ELSE
	142	CALL iom_put( "ssh" , sshn ) ! sea surface height
	143	ENDIF
	144
[7646]	145	IF( iom_use("wetdep") ) & ! wet depth
[9023]	146	CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
[5107]	147
	148	CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature
	149	CALL iom_put( "sst", tsn(:,:,1,jp_tem) ) ! surface temperature
	150	IF ( iom_use("sbt") ) THEN
[4990]	151	DO jj = 1, jpj
	152	DO ji = 1, jpi
[9019]	153	ikbot = mbkt(ji,jj)
	154	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
[4990]	155	END DO
[5107]	156	END DO
	157	CALL iom_put( "sbt", z2d ) ! bottom temperature
	158	ENDIF
	159
	160	CALL iom_put( "soce", tsn(:,:,:,jp_sal) ) ! 3D salinity
	161	CALL iom_put( "sss", tsn(:,:,1,jp_sal) ) ! surface salinity
	162	IF ( iom_use("sbs") ) THEN
[4990]	163	DO jj = 1, jpj
	164	DO ji = 1, jpi
[9019]	165	ikbot = mbkt(ji,jj)
	166	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
[4990]	167	END DO
[5107]	168	END DO
	169	CALL iom_put( "sbs", z2d ) ! bottom salinity
	170	ENDIF
[5463]	171
	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	!
[6140]	394
[9019]	395	IF (ln_diatmb) CALL dia_tmb ! tmb values
	396
	397	IF (ln_dia25h) CALL dia_25h( kt ) ! 25h averaging
[6140]	398
[9124]	399	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	400	!
[1482]	401	END SUBROUTINE dia_wri
	402
	403	#else
[2528]	404	!!----------------------------------------------------------------------
	405	!! Default option use IOIPSL library
	406	!!----------------------------------------------------------------------
	407
[9652]	408	INTEGER FUNCTION dia_wri_alloc()
	409	!!----------------------------------------------------------------------
	410	INTEGER, DIMENSION(2) :: ierr
	411	!!----------------------------------------------------------------------
	412	ierr = 0
	413	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	414	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	415	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
	416	!
	417	dia_wri_alloc = MAXVAL(ierr)
[10425]	418	CALL mpp_sum( 'diawri', dia_wri_alloc )
[9652]	419	!
	420	END FUNCTION dia_wri_alloc
	421
	422
[1561]	423	SUBROUTINE dia_wri( kt )
[3]	424	!!---------------------------------------------------------------------
	425	!! * ROUTINE dia_wri *
	426	!!
	427	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	428	!! NETCDF format is used by default
	429	!!
	430	!! ** Method : At the beginning of the first time step (nit000),
	431	!! define all the NETCDF files and fields
	432	!! At each time step call histdef to compute the mean if ncessary
[11536]	433	!! Each nn_write time step, output the instantaneous or mean fields
[3]	434	!!----------------------------------------------------------------------
[5836]	435	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	436	!
[2528]	437	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	438	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	439	INTEGER :: inum = 11 ! temporary logical unit
[3294]	440	INTEGER :: ji, jj, jk ! dummy loop indices
	441	INTEGER :: ierr ! error code return from allocation
[2528]	442	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[3609]	443	INTEGER :: jn, ierror ! local integers
[6140]	444	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	445	!
[9019]	446	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
	447	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
[3]	448	!!----------------------------------------------------------------------
[1482]	449	!
[9019]	450	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[10425]	451	CALL dia_wri_state( 'output.init' )
[1482]	452	ninist = 0
	453	ENDIF
	454	!
[11536]	455	IF( nn_write == -1 ) RETURN ! we will never do any output
	456	!
	457	IF( ln_timing ) CALL timing_start('dia_wri')
	458	!
[3]	459	! 0. Initialisation
	460	! -----------------
[632]	461
[9019]	462	ll_print = .FALSE. ! local variable for debugging
[3]	463	ll_print = ll_print .AND. lwp
	464
	465	! Define frequency of output and means
[5566]	466	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]	467	#if defined key_diainstant
[11536]	468	zsto = nn_write * rdt
[1312]	469	clop = "inst("//TRIM(clop)//")"
[3]	470	#else
[6140]	471	zsto=rdt
[1312]	472	clop = "ave("//TRIM(clop)//")"
[3]	473	#endif
[11536]	474	zout = nn_write * rdt
[6140]	475	zmax = ( nitend - nit000 + 1 ) * rdt
[3]	476
	477	! Define indices of the horizontal output zoom and vertical limit storage
	478	iimi = 1 ; iima = jpi
	479	ijmi = 1 ; ijma = jpj
	480	ipk = jpk
	481
	482	! define time axis
[1334]	483	it = kt
	484	itmod = kt - nit000 + 1
[3]	485
	486
	487	! 1. Define NETCDF files and fields at beginning of first time step
	488	! -----------------------------------------------------------------
	489
	490	IF( kt == nit000 ) THEN
	491
	492	! Define the NETCDF files (one per grid)
[632]	493
[3]	494	! Compute julian date from starting date of the run
[1309]	495	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	496	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	497	IF(lwp)WRITE(numout,*)
	498	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	499	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	500	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	501	' limit storage in depth = ', ipk
	502
	503	! WRITE root name in date.file for use by postpro
[1581]	504	IF(lwp) THEN
[11536]	505	CALL dia_nam( clhstnam, nn_write,' ' )
[1581]	506	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	507	WRITE(inum,*) clhstnam
	508	CLOSE(inum)
	509	ENDIF
[632]	510
[3]	511	! Define the T grid FILE ( nid_T )
[632]	512
[11536]	513	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[3]	514	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	515	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	516	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	517	& nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	518	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	519	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	520	! ! Index of ocean points
	521	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	522	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	523	!
	524	IF( ln_icebergs ) THEN
	525	!
	526	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	527	!! that routine is called from nemogcm, so do it here immediately before its needed
	528	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	529	CALL mpp_sum( 'diawri', ierror )
[3609]	530	IF( ierror /= 0 ) THEN
	531	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	532	RETURN
	533	ENDIF
	534	!
	535	!! iceberg vertical coordinate is class number
	536	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	537	& "number", nclasses, class_num, nb_T )
	538	!
	539	!! each class just needs the surface index pattern
	540	ndim_bT = 3
	541	DO jn = 1,nclasses
	542	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	543	ENDDO
	544	!
	545	ENDIF
[3]	546
	547	! Define the U grid FILE ( nid_U )
	548
[11536]	549	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[3]	550	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	551	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	552	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	553	& nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	554	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	555	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	556	! ! Index of ocean points
	557	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	558	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	559
	560	! Define the V grid FILE ( nid_V )
	561
[11536]	562	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[3]	563	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	564	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	565	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	566	& nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	567	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	568	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	569	! ! Index of ocean points
	570	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	571	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	572
	573	! Define the W grid FILE ( nid_W )
	574
[11536]	575	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[3]	576	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	577	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	578	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	579	& nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	580	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	581	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	582
[632]	583
[3]	584	! Declare all the output fields as NETCDF variables
	585
	586	! !!! nid_T : 3D
	587	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	588	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	589	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	590	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]	591	IF( .NOT.ln_linssh ) THEN
[4292]	592	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	593	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	594	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	595	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	596	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	597	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	598	ENDIF
[3]	599	! !!! nid_T : 2D
	600	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	601	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	602	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	603	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	604	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	605	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	606	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	607	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	608	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	609	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	610	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	611	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[6140]	612	IF( ln_linssh ) THEN
[4292]	613	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
[3625]	614	& , "KgC/m2/s", & ! sosst_cd
[4292]	615	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	616	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
[3625]	617	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	618	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	619	ENDIF
[888]	620	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	621	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	622	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	623	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	624	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	625	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	626	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	627	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	628	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	629	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	630	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	631	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3609]	632	!
	633	IF( ln_icebergs ) THEN
	634	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	635	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	636	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	637	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	638	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	639	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	640	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	641	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	642	IF( ln_bergdia ) THEN
	643	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	644	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	645	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	646	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	647	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	648	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	649	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	650	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	651	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	652	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	653	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	654	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	655	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	656	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	657	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	658	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	659	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	660	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	661	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	662	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	663	ENDIF
	664	ENDIF
	665
[11536]	666	IF( ln_ssr ) THEN
[4990]	667	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	668	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	669	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	670	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	671	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	672	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	673	ENDIF
[11536]	674
[3]	675	clmx ="l_max(only(x))" ! max index on a period
[5836]	676	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	677	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
[3]	678	#if defined key_diahth
	679	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	680	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	681	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	682	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	683	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	684	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[7646]	685	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
[3]	686	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	687	#endif
	688
[2528]	689	CALL histend( nid_T, snc4chunks=snc4set )
[3]	690
	691	! !!! nid_U : 3D
	692	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	693	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]	694	IF( ln_wave .AND. ln_sdw) THEN
	695	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	696	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	697	ENDIF
[3]	698	! !!! nid_U : 2D
[888]	699	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	700	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	701
[2528]	702	CALL histend( nid_U, snc4chunks=snc4set )
[3]	703
	704	! !!! nid_V : 3D
	705	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	706	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]	707	IF( ln_wave .AND. ln_sdw) THEN
	708	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	709	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	710	ENDIF
[3]	711	! !!! nid_V : 2D
[888]	712	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	713	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	714
[2528]	715	CALL histend( nid_V, snc4chunks=snc4set )
[3]	716
	717	! !!! nid_W : 3D
	718	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	719	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	720	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	721	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]	722	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
[255]	723	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	724
[9019]	725	IF( ln_zdfddm ) THEN
[3]	726	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	727	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	728	ENDIF
[7646]	729
	730	IF( ln_wave .AND. ln_sdw) THEN
	731	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	732	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	733	ENDIF
[3]	734	! !!! nid_W : 2D
[2528]	735	CALL histend( nid_W, snc4chunks=snc4set )
[3]	736
	737	IF(lwp) WRITE(numout,*)
	738	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	739	IF(ll_print) CALL FLUSH(numout )
	740
	741	ENDIF
	742
	743	! 2. Start writing data
	744	! ---------------------
	745
[4292]	746	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	747	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	748	! donne le nombre d'elements, et ndex la liste des indices a sortir
	749
[11536]	750	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
[3]	751	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	752	WRITE(numout,*) '~~~~~~ '
	753	ENDIF
	754
[6140]	755	IF( .NOT.ln_linssh ) THEN
	756	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	757	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	758	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	759	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	760	ELSE
	761	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	762	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	763	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	764	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	765	ENDIF
[6140]	766	IF( .NOT.ln_linssh ) THEN
[7753]	767	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
[6140]	768	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	769	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
[4292]	770	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	771	ENDIF
[3]	772	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
[2528]	773	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	774	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	775	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	776	! (includes virtual salt flux beneath ice
	777	! in linear free surface case)
[6140]	778	IF( ln_linssh ) THEN
[7753]	779	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
[4292]	780	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
[7753]	781	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
[4292]	782	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	783	ENDIF
[888]	784	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	785	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	786	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	787	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	788	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	789	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[3609]	790	!
	791	IF( ln_icebergs ) THEN
	792	!
	793	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	794	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	795	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	796	!
	797	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	798	!
	799	IF( ln_bergdia ) THEN
	800	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	801	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	802	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	803	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	804	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	805	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	806	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	807	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	808	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	809	!
	810	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	811	ENDIF
	812	ENDIF
	813
[11536]	814	IF( ln_ssr ) THEN
[4990]	815	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	816	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[11536]	817	zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	818	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	819	ENDIF
	820	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	821	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	822
	823	#if defined key_diahth
	824	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	825	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	826	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	827	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	828	#endif
[888]	829
[3]	830	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
[888]	831	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	832
	833	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
[888]	834	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	835
[11418]	836	IF( ln_zad_Aimp ) THEN
	837	CALL histwrite( nid_W, "vovecrtz", it, wn + wi , ndim_T, ndex_T ) ! vert. current
	838	ELSE
	839	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
	840	ENDIF
[3]	841	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	842	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	843	IF( ln_zdfddm ) THEN
	844	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	845	ENDIF
	846
[7646]	847	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	848	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	849	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	850	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	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	!
[9124]	862	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	863	!
[3]	864	END SUBROUTINE dia_wri
[1567]	865	#endif
	866
[10425]	867	SUBROUTINE dia_wri_state( cdfile_name )
[3]	868	!!---------------------------------------------------------------------
	869	!! * ROUTINE dia_wri_state *
	870	!!
	871	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	872	!! the instantaneous ocean state and forcing fields.
	873	!! Used to find errors in the initial state or save the last
	874	!! ocean state in case of abnormal end of a simulation
	875	!!
	876	!! ** Method : NetCDF files using ioipsl
	877	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	878	!! File 'output.abort.nc' is created in case of abnormal job end
	879	!!----------------------------------------------------------------------
[1334]	880	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	881	!!
	882	INTEGER :: inum
[3]	883	!!----------------------------------------------------------------------
[3294]	884	!
[648]	885	IF(lwp) WRITE(numout,*)
	886	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	887	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
[10425]	888	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
[648]	889
[9570]	890	#if defined key_si3
[10425]	891	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
[1482]	892	#else
[10425]	893	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
[1482]	894	#endif
[3]	895
[10425]	896	CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) ) ! now temperature
	897	CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) ) ! now salinity
	898	CALL iom_rstput( 0, 0, inum, 'sossheig', sshn ) ! sea surface height
	899	CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity
	900	CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity
[11418]	901	IF( ln_zad_Aimp ) THEN
	902	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi ) ! now k-velocity
	903	ELSE
	904	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity
	905	ENDIF
[7646]	906	IF( ALLOCATED(ahtu) ) THEN
[10425]	907	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	908	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	909	ENDIF
	910	IF( ALLOCATED(ahmt) ) THEN
[10425]	911	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	912	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	913	ENDIF
[10425]	914	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	915	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	916	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	917	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	918	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	919	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[6140]	920	IF( .NOT.ln_linssh ) THEN
[10425]	921	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
	922	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
	923	END IF
[7646]	924	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	925	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	926	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	927	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	928	ENDIF
[12386]	929
	930	IF( ln_zdfosm ) THEN
	931	CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1) ) ! now boundary-layer depth
	932	CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1) ) ! now mixed-layer depth
	933	CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask ) ! w-level diffusion
	934	CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask ) ! now w-level viscosity
	935	CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask ) ! non-local t forcing
	936	CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask ) ! non-local s forcing
	937	CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*wmask ) ! non-local u forcing
	938	CALL iom_rstput( 0, 0, inum, 'ghamv', ghamu*wmask ) ! non-local v forcing
	939	IF( ln_osm_mle ) THEN
	940	CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1) ) ! now transition-layer depth
	941	END IF
	942	ENDIF
[10425]	943
	944	#if defined key_si3
	945	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
	946	CALL ice_wri_state( inum )
[1561]	947	ENDIF
	948	#endif
[10425]	949	!
	950	CALL iom_close( inum )
[3294]	951	!
[3]	952	END SUBROUTINE dia_wri_state
[6140]	953
[3]	954	!!======================================================================
	955	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/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/DIA/diawri.F90 @ 12386

Download in other formats: