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/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90 @ 12667

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