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/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90 @ 11348

Last change on this file since 11348 was 11348, checked in by gsamson, 5 years ago

dev_r11265_ABL :

merge HPC-13_IRRMANN_BDY_optimization branch @ rev11332 with dev_r11265_ABL branch @ rev11334
allow ln_dm2dc option with ABL
cosmetic change in sbcabl.F90

identical results with rev11334 for bulk and abl orca2

Property svn:keywords set to Id

File size: 53.5 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
[11334]	28	USE abl ! abl variables in case ln_abl = .true.
[9019]	29	USE dom_oce ! ocean space and time domain
	30	USE phycst ! physical constants
	31	USE dianam ! build name of file (routine)
	32	USE diahth ! thermocline diagnostics
	33	USE dynadv , ONLY: ln_dynadv_vec
	34	USE icb_oce ! Icebergs
	35	USE icbdia ! Iceberg budgets
	36	USE ldftra ! lateral physics: eddy diffusivity coef.
	37	USE ldfdyn ! lateral physics: eddy viscosity coef.
	38	USE sbc_oce ! Surface boundary condition: ocean fields
	39	USE sbc_ice ! Surface boundary condition: ice fields
	40	USE sbcssr ! restoring term toward SST/SSS climatology
	41	USE sbcwave ! wave parameters
	42	USE wet_dry ! wetting and drying
	43	USE zdf_oce ! ocean vertical physics
	44	USE zdfdrg ! ocean vertical physics: top/bottom friction
	45	USE zdfmxl ! 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
[11334]	69	#if ! defined key_iomput
	70	PUBLIC dia_wri_alloc_abl ! Called by sbcabl module (if ln_abl = .true.)
	71	#endif
[2528]	72	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[11334]	73	INTEGER :: nb_T, ndim_bT ! grid_T file
[2528]	74	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	75	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	76	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
[11334]	77	INTEGER :: ndim_A, ndim_hA ! grid_T file
	78	INTEGER :: nid_A, nz_A, nh_A ! grid_ABL file
[2528]	79	INTEGER :: ndex(1) ! ???
[11334]	80	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV, ndex_hA
	81	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V, ndex_A
[3609]	82	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	83
	84	!! * Substitutions
[1756]	85	# include "vectopt_loop_substitute.h90"
[3]	86	!!----------------------------------------------------------------------
[9598]	87	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	88	!! $Id$
[10068]	89	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	90	!!----------------------------------------------------------------------
	91	CONTAINS
	92
[6140]	93	#if defined key_iomput
[3]	94	!!----------------------------------------------------------------------
[2528]	95	!! 'key_iomput' use IOM library
	96	!!----------------------------------------------------------------------
[9652]	97	INTEGER FUNCTION dia_wri_alloc()
	98	!
	99	dia_wri_alloc = 0
	100	!
	101	END FUNCTION dia_wri_alloc
[2715]	102
[9652]	103
[1561]	104	SUBROUTINE dia_wri( kt )
[1482]	105	!!---------------------------------------------------------------------
	106	!! * ROUTINE dia_wri *
	107	!!
	108	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	109	!! NETCDF format is used by default
	110	!!
	111	!! ** Method : use iom_put
	112	!!----------------------------------------------------------------------
	113	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1756]	114	!!
[9019]	115	INTEGER :: ji, jj, jk ! dummy loop indices
	116	INTEGER :: ikbot ! local integer
	117	REAL(wp):: zztmp , zztmpx ! local scalar
	118	REAL(wp):: zztmp2, zztmpy ! - -
	119	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
	120	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
[1482]	121	!!----------------------------------------------------------------------
	122	!
[9124]	123	IF( ln_timing ) CALL timing_start('dia_wri')
[3294]	124	!
[1482]	125	! Output the initial state and forcings
	126	IF( ninist == 1 ) THEN
[10425]	127	CALL dia_wri_state( 'output.init' )
[1482]	128	ninist = 0
	129	ENDIF
[3]	130
[6351]	131	! Output of initial vertical scale factor
	132	CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10114]	133	CALL iom_put("e3u_0", e3u_0(:,:,:) )
	134	CALL iom_put("e3v_0", e3v_0(:,:,:) )
[6351]	135	!
[6387]	136	CALL iom_put( "e3t" , e3t_n(:,:,:) )
	137	CALL iom_put( "e3u" , e3u_n(:,:,:) )
	138	CALL iom_put( "e3v" , e3v_n(:,:,:) )
	139	CALL iom_put( "e3w" , e3w_n(:,:,:) )
[6351]	140	IF( iom_use("e3tdef") ) &
[6387]	141	CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
[5461]	142
[9023]	143	IF( ll_wd ) THEN
	144	CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)
	145	ELSE
	146	CALL iom_put( "ssh" , sshn ) ! sea surface height
	147	ENDIF
	148
[7646]	149	IF( iom_use("wetdep") ) & ! wet depth
[9023]	150	CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
[5107]	151
	152	CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature
	153	CALL iom_put( "sst", tsn(:,:,1,jp_tem) ) ! surface temperature
	154	IF ( iom_use("sbt") ) THEN
[4990]	155	DO jj = 1, jpj
	156	DO ji = 1, jpi
[9019]	157	ikbot = mbkt(ji,jj)
	158	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
[4990]	159	END DO
[5107]	160	END DO
	161	CALL iom_put( "sbt", z2d ) ! bottom temperature
	162	ENDIF
	163
	164	CALL iom_put( "soce", tsn(:,:,:,jp_sal) ) ! 3D salinity
	165	CALL iom_put( "sss", tsn(:,:,1,jp_sal) ) ! surface salinity
	166	IF ( iom_use("sbs") ) THEN
[4990]	167	DO jj = 1, jpj
	168	DO ji = 1, jpi
[9019]	169	ikbot = mbkt(ji,jj)
	170	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
[4990]	171	END DO
[5107]	172	END DO
	173	CALL iom_put( "sbs", z2d ) ! bottom salinity
	174	ENDIF
[5463]	175
	176	IF ( iom_use("taubot") ) THEN ! bottom stress
[9019]	177	zztmp = rau0 * 0.25
[7753]	178	z2d(:,:) = 0._wp
[5463]	179	DO jj = 2, jpjm1
	180	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	181	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * un(ji ,jj,mbku(ji ,jj)) )**2 &
	182	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj)) )**2 &
	183	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vn(ji,jj ,mbkv(ji,jj )) )**2 &
	184	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1)) )**2
	185	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
[5463]	186	!
[9019]	187	END DO
	188	END DO
[10425]	189	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
[5463]	190	CALL iom_put( "taubot", z2d )
	191	ENDIF
[5107]	192
[9019]	193	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
	194	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
[5107]	195	IF ( iom_use("sbu") ) THEN
[4990]	196	DO jj = 1, jpj
	197	DO ji = 1, jpi
[9019]	198	ikbot = mbku(ji,jj)
	199	z2d(ji,jj) = un(ji,jj,ikbot)
[4990]	200	END DO
[5107]	201	END DO
	202	CALL iom_put( "sbu", z2d ) ! bottom i-current
	203	ENDIF
	204
[9019]	205	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
	206	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
[5107]	207	IF ( iom_use("sbv") ) THEN
[4990]	208	DO jj = 1, jpj
	209	DO ji = 1, jpi
[9019]	210	ikbot = mbkv(ji,jj)
	211	z2d(ji,jj) = vn(ji,jj,ikbot)
[4990]	212	END DO
[5107]	213	END DO
	214	CALL iom_put( "sbv", z2d ) ! bottom j-current
[4990]	215	ENDIF
[1482]	216
[5461]	217	CALL iom_put( "woce", wn ) ! vertical velocity
	218	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
	219	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[7753]	220	z2d(:,:) = rau0 * e1e2t(:,:)
[5461]	221	DO jk = 1, jpk
[7753]	222	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
[5461]	223	END DO
	224	CALL iom_put( "w_masstr" , z3d )
	225	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
	226	ENDIF
	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) )
[11334]	416	dia_wri_alloc = MAXVAL(ierr)
[10425]	417	CALL mpp_sum( 'diawri', dia_wri_alloc )
[9652]	418	!
	419	END FUNCTION dia_wri_alloc
[11334]	420
	421	INTEGER FUNCTION dia_wri_alloc_abl()
	422	!!----------------------------------------------------------------------
	423	ALLOCATE( ndex_hA(jpijpj), ndex_A (jpijpj*jpkam1), STAT=dia_wri_alloc_abl)
	424	CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
	425	!
	426	END FUNCTION dia_wri_alloc_abl
	427
[1561]	428	SUBROUTINE dia_wri( kt )
[3]	429	!!---------------------------------------------------------------------
	430	!! * ROUTINE dia_wri *
	431	!!
	432	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	433	!! NETCDF format is used by default
	434	!!
	435	!! ** Method : At the beginning of the first time step (nit000),
	436	!! define all the NETCDF files and fields
	437	!! At each time step call histdef to compute the mean if ncessary
	438	!! Each nwrite time step, output the instantaneous or mean fields
	439	!!----------------------------------------------------------------------
[5836]	440	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	441	!
[2528]	442	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	443	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	444	INTEGER :: inum = 11 ! temporary logical unit
[3294]	445	INTEGER :: ji, jj, jk ! dummy loop indices
	446	INTEGER :: ierr ! error code return from allocation
[11334]	447	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma, ipka ! local integers
[3609]	448	INTEGER :: jn, ierror ! local integers
[6140]	449	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	450	!
[9019]	451	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
	452	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
[11334]	453	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! 3D workspace
[3]	454	!!----------------------------------------------------------------------
[3294]	455	!
[9124]	456	IF( ln_timing ) CALL timing_start('dia_wri')
[1482]	457	!
[9019]	458	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[10425]	459	CALL dia_wri_state( 'output.init' )
[1482]	460	ninist = 0
	461	ENDIF
[11334]	462
[1482]	463	!
[3]	464	! 0. Initialisation
	465	! -----------------
[632]	466
[9019]	467	ll_print = .FALSE. ! local variable for debugging
[3]	468	ll_print = ll_print .AND. lwp
	469
	470	! Define frequency of output and means
[5566]	471	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]	472	#if defined key_diainstant
[6140]	473	zsto = nwrite * rdt
[1312]	474	clop = "inst("//TRIM(clop)//")"
[3]	475	#else
[6140]	476	zsto=rdt
[1312]	477	clop = "ave("//TRIM(clop)//")"
[3]	478	#endif
[6140]	479	zout = nwrite * rdt
	480	zmax = ( nitend - nit000 + 1 ) * rdt
[3]	481
	482	! Define indices of the horizontal output zoom and vertical limit storage
	483	iimi = 1 ; iima = jpi
	484	ijmi = 1 ; ijma = jpj
	485	ipk = jpk
[11334]	486	IF(ln_abl) ipka = jpkam1
[3]	487
	488	! define time axis
[1334]	489	it = kt
	490	itmod = kt - nit000 + 1
[3]	491
	492
	493	! 1. Define NETCDF files and fields at beginning of first time step
	494	! -----------------------------------------------------------------
	495
	496	IF( kt == nit000 ) THEN
	497
	498	! Define the NETCDF files (one per grid)
[632]	499
[3]	500	! Compute julian date from starting date of the run
[1309]	501	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	502	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	503	IF(lwp)WRITE(numout,*)
	504	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	505	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	506	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	507	' limit storage in depth = ', ipk
	508
	509	! WRITE root name in date.file for use by postpro
[1581]	510	IF(lwp) THEN
[895]	511	CALL dia_nam( clhstnam, nwrite,' ' )
[1581]	512	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	513	WRITE(inum,*) clhstnam
	514	CLOSE(inum)
	515	ENDIF
[632]	516
[3]	517	! Define the T grid FILE ( nid_T )
[632]	518
[3]	519	CALL dia_nam( clhstnam, nwrite, 'grid_T' )
	520	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	521	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	522	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	523	& nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	524	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	525	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	526	! ! Index of ocean points
	527	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	528	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	529	!
	530	IF( ln_icebergs ) THEN
	531	!
	532	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	533	!! that routine is called from nemogcm, so do it here immediately before its needed
	534	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	535	CALL mpp_sum( 'diawri', ierror )
[3609]	536	IF( ierror /= 0 ) THEN
	537	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	538	RETURN
	539	ENDIF
	540	!
	541	!! iceberg vertical coordinate is class number
	542	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	543	& "number", nclasses, class_num, nb_T )
	544	!
	545	!! each class just needs the surface index pattern
	546	ndim_bT = 3
	547	DO jn = 1,nclasses
	548	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	549	ENDDO
	550	!
	551	ENDIF
[3]	552
	553	! Define the U grid FILE ( nid_U )
	554
	555	CALL dia_nam( clhstnam, nwrite, 'grid_U' )
	556	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	557	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	558	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	559	& nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	560	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	561	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	562	! ! Index of ocean points
	563	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	564	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	565
	566	! Define the V grid FILE ( nid_V )
	567
	568	CALL dia_nam( clhstnam, nwrite, 'grid_V' ) ! filename
	569	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	570	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	571	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	572	& nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	573	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	574	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	575	! ! Index of ocean points
	576	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	577	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	578
	579	! Define the W grid FILE ( nid_W )
	580
	581	CALL dia_nam( clhstnam, nwrite, 'grid_W' ) ! filename
	582	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	583	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	584	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[6140]	585	& nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	586	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	587	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	588
[11334]	589	IF( ln_abl ) THEN
	590	! Define the ABL grid FILE ( nid_A )
	591	CALL dia_nam( clhstnam, nwrite, 'grid_ABL' )
	592	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	593	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	594	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
	595	& nit000-1, zjulian, rdt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
	596	CALL histvert( nid_A, "ght_abl", "Vertical T levels", & ! Vertical grid: gdept
	597	& "m", ipka, ght_abl(2:jpka), nz_A, "up" )
	598	! ! Index of ocean points
	599	ALLOCATE( zw3d_abl(jpi,jpj,ipka) )
	600	zw3d_abl(:,:,:) = 1._wp
	601	CALL wheneq( jpijpjipka, zw3d_abl, 1, 1., ndex_A , ndim_A ) ! volume
	602	CALL wheneq( jpi*jpj , zw3d_abl, 1, 1., ndex_hA, ndim_hA ) ! surface
	603	DEALLOCATE(zw3d_abl)
	604	ENDIF
[632]	605
[3]	606	! Declare all the output fields as NETCDF variables
	607
	608	! !!! nid_T : 3D
	609	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	610	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	611	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	612	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]	613	IF( .NOT.ln_linssh ) THEN
[4292]	614	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	615	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	616	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	617	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	618	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	619	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	620	ENDIF
[3]	621	! !!! nid_T : 2D
	622	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	623	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	624	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	625	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	626	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	627	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	628	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	629	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	630	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	631	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	632	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	633	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[6140]	634	IF( ln_linssh ) THEN
[4292]	635	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
[3625]	636	& , "KgC/m2/s", & ! sosst_cd
[4292]	637	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	638	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
[3625]	639	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	640	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	641	ENDIF
[888]	642	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	643	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	644	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	645	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	646	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	647	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	648	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	649	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	650	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	651	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	652	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	653	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3609]	654	!
[11334]	655	IF( ln_abl ) THEN
	656	! !!! nid_A : 3D
	657	CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl
	658	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	659	CALL histdef( nid_A, "q_abl", "Humidity" , "kg/kg" , & ! q_abl
	660	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	661	CALL histdef( nid_A, "u_abl", "Atmospheric U-wind " , "m/s" , & ! u_abl
	662	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	663	CALL histdef( nid_A, "v_abl", "Atmospheric V-wind " , "m/s" , & ! v_abl
	664	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	665	CALL histdef( nid_A, "tke_abl", "Atmospheric TKE " , "m2/s2" , & ! tke_abl
	666	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	667	CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s" , & ! avm_abl
	668	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	669	CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2", & ! avt_abl
	670	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	671	CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height " , "m", & ! pblh
	672	& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	673	#if defined key_si3
	674	CALL histdef( nid_A, "oce_frac", "Fraction of open ocean" , " ", & ! ato_i
	675	& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	676	#endif
	677	CALL histend( nid_A, snc4chunks=snc4set )
	678	!
	679	ENDIF
	680	!
[3609]	681	IF( ln_icebergs ) THEN
	682	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	683	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	684	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	685	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	686	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	687	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	688	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	689	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	690	IF( ln_bergdia ) THEN
	691	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	692	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	693	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	694	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	695	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	696	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	697	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	698	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	699	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	700	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	701	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	702	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	703	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	704	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	705	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	706	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	707	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	708	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	709	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	710	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	711	ENDIF
	712	ENDIF
	713
[11348]	714	IF( ln_ssr ) THEN
[4990]	715	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	716	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	717	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	718	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	719	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	720	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	721	ENDIF
[11348]	722
[3]	723	clmx ="l_max(only(x))" ! max index on a period
[11334]	724	CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	725	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
[3]	726	#if defined key_diahth
	727	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	728	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	729	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	730	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	731	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	732	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[7646]	733	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
[3]	734	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	735	#endif
	736
[2528]	737	CALL histend( nid_T, snc4chunks=snc4set )
[3]	738
	739	! !!! nid_U : 3D
	740	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	741	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]	742	IF( ln_wave .AND. ln_sdw) THEN
	743	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	744	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	745	ENDIF
[3]	746	! !!! nid_U : 2D
[888]	747	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	748	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	749
[2528]	750	CALL histend( nid_U, snc4chunks=snc4set )
[3]	751
	752	! !!! nid_V : 3D
	753	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	754	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]	755	IF( ln_wave .AND. ln_sdw) THEN
	756	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	757	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	758	ENDIF
[3]	759	! !!! nid_V : 2D
[888]	760	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	761	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	762
[2528]	763	CALL histend( nid_V, snc4chunks=snc4set )
[3]	764
	765	! !!! nid_W : 3D
	766	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	767	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	768	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	769	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]	770	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
[255]	771	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	772
[9019]	773	IF( ln_zdfddm ) THEN
[3]	774	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	775	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	776	ENDIF
[7646]	777
	778	IF( ln_wave .AND. ln_sdw) THEN
	779	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	780	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	781	ENDIF
[3]	782	! !!! nid_W : 2D
[2528]	783	CALL histend( nid_W, snc4chunks=snc4set )
[3]	784
	785	IF(lwp) WRITE(numout,*)
	786	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	787	IF(ll_print) CALL FLUSH(numout )
	788
	789	ENDIF
	790
	791	! 2. Start writing data
	792	! ---------------------
	793
[4292]	794	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	795	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	796	! donne le nombre d'elements, et ndex la liste des indices a sortir
	797
[1334]	798	IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
[3]	799	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	800	WRITE(numout,*) '~~~~~~ '
	801	ENDIF
	802
[6140]	803	IF( .NOT.ln_linssh ) THEN
	804	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	805	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	806	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	807	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	808	ELSE
	809	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	810	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	811	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	812	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	813	ENDIF
[6140]	814	IF( .NOT.ln_linssh ) THEN
[7753]	815	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
[6140]	816	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	817	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
[4292]	818	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	819	ENDIF
[3]	820	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
[2528]	821	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	822	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	823	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	824	! (includes virtual salt flux beneath ice
	825	! in linear free surface case)
[6140]	826	IF( ln_linssh ) THEN
[7753]	827	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
[4292]	828	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
[7753]	829	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
[4292]	830	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	831	ENDIF
[888]	832	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	833	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	834	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	835	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	836	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	837	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[3609]	838	!
[11334]	839	IF( ln_abl ) THEN
	840	ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
	841	IF( ln_mskland ) THEN
	842	DO jk=1,jpka
	843	zw3d_abl(:,:,jk) = tmask(:,:,1)
	844	END DO
	845	ELSE
	846	zw3d_abl(:,:,:) = 1._wp
	847	ENDIF
	848	CALL histwrite( nid_A, "pblh" , it, pblh(:,:) *zw3d_abl(:,:,1 ), ndim_hA, ndex_hA ) ! pblh
	849	CALL histwrite( nid_A, "u_abl" , it, u_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! u_abl
	850	CALL histwrite( nid_A, "v_abl" , it, v_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! v_abl
	851	CALL histwrite( nid_A, "t_abl" , it, tq_abl (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! t_abl
	852	CALL histwrite( nid_A, "q_abl" , it, tq_abl (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! q_abl
	853	CALL histwrite( nid_A, "tke_abl", it, tke_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! tke_abl
	854	CALL histwrite( nid_A, "avm_abl", it, avm_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avm_abl
	855	CALL histwrite( nid_A, "avt_abl", it, avt_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avt_abl
	856	#if defined key_si3
	857	CALL histwrite( nid_A, "oce_frac" , it, ato_i(:,:) , ndim_hA, ndex_hA ) ! ato_i
	858	#endif
	859	DEALLOCATE(zw3d_abl)
	860	ENDIF
	861	!
[3609]	862	IF( ln_icebergs ) THEN
	863	!
	864	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	865	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	866	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	867	!
	868	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	869	!
	870	IF( ln_bergdia ) THEN
	871	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	872	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	873	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	874	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	875	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	876	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	877	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	878	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	879	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	880	!
	881	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	882	ENDIF
	883	ENDIF
	884
[11348]	885	IF( ln_ssr ) THEN
[4990]	886	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	887	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[11348]	888	zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	889	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	890	ENDIF
[11334]	891	zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	892	CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	893
	894	#if defined key_diahth
	895	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	896	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	897	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	898	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	899	#endif
[888]	900
[3]	901	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
[888]	902	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	903
	904	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
[888]	905	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	906
	907	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
	908	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	909	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	910	IF( ln_zdfddm ) THEN
	911	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	912	ENDIF
	913
[7646]	914	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	915	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	916	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	917	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	918	ENDIF
	919
[1318]	920	! 3. Close all files
	921	! ---------------------------------------
[1561]	922	IF( kt == nitend ) THEN
[3]	923	CALL histclo( nid_T )
	924	CALL histclo( nid_U )
	925	CALL histclo( nid_V )
	926	CALL histclo( nid_W )
[11334]	927	IF(ln_abl) CALL histclo( nid_A )
[3]	928	ENDIF
[2528]	929	!
[9124]	930	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	931	!
[3]	932	END SUBROUTINE dia_wri
[1567]	933	#endif
	934
[10425]	935	SUBROUTINE dia_wri_state( cdfile_name )
[3]	936	!!---------------------------------------------------------------------
	937	!! * ROUTINE dia_wri_state *
	938	!!
	939	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	940	!! the instantaneous ocean state and forcing fields.
	941	!! Used to find errors in the initial state or save the last
	942	!! ocean state in case of abnormal end of a simulation
	943	!!
	944	!! ** Method : NetCDF files using ioipsl
	945	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	946	!! File 'output.abort.nc' is created in case of abnormal job end
	947	!!----------------------------------------------------------------------
[1334]	948	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	949	!!
	950	INTEGER :: inum
[3]	951	!!----------------------------------------------------------------------
[3294]	952	!
[648]	953	IF(lwp) WRITE(numout,*)
	954	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	955	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
[10425]	956	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
[648]	957
[9570]	958	#if defined key_si3
[10425]	959	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
[1482]	960	#else
[10425]	961	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
[1482]	962	#endif
[3]	963
[10425]	964	CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) ) ! now temperature
	965	CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) ) ! now salinity
	966	CALL iom_rstput( 0, 0, inum, 'sossheig', sshn ) ! sea surface height
	967	CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity
	968	CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity
	969	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity
[7646]	970	IF( ALLOCATED(ahtu) ) THEN
[10425]	971	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	972	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	973	ENDIF
	974	IF( ALLOCATED(ahmt) ) THEN
[10425]	975	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	976	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	977	ENDIF
[10425]	978	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	979	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	980	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	981	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	982	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	983	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[6140]	984	IF( .NOT.ln_linssh ) THEN
[10425]	985	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
	986	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
	987	END IF
[7646]	988	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	989	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	990	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	991	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	992	ENDIF
[10425]	993
	994	#if defined key_si3
	995	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
	996	CALL ice_wri_state( inum )
[1561]	997	ENDIF
	998	#endif
[10425]	999	!
	1000	CALL iom_close( inum )
[3294]	1001	!
[3]	1002	END SUBROUTINE dia_wri_state
[6140]	1003
[3]	1004	!!======================================================================
	1005	END MODULE diawri

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

Download in other formats: