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/UKMO/NEMO_4.0.3_GM_rossby_radius_ramp/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.3_GM_rossby_radius_ramp/src/OCE/DIA/diawri.F90 @ 13883

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

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

Download in other formats: