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 branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6005

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

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

Download in other formats: