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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90 @ 11822

Last change on this file since 11822 was 11822, checked in by acc, 4 years ago
Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.
Property svn:keywords set to `Id`
File size: 49.7 KB

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

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

Download in other formats: