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_r11943_MERGE_2019/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diawri.F90 @ 12150

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

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

Download in other formats: