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 trunk/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6667

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

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

Download in other formats: