Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

diaopfoam.F90

Last change on this file was 12918, checked in by petesykes, 4 years ago
Land masking for operational diagnostics
File size: 19.9 KB

Line
1	MODULE diaopfoam
2	!!======================================================================
3	!! * MODULE diaopfoam *
4	!! Output stream for operational use
5	!!======================================================================
6	!! History : 3.6 ! 2016 (P Sykes) Original code
7	!!----------------------------------------------------------------------
8	USE oce ! ocean dynamics and tracers variables
9	USE dom_oce ! ocean space and time domain
10	USE diainsitutem, ONLY: rinsitu_t, theta2t
11	USE in_out_manager ! I/O units
12	USE iom ! I/0 library
13	USE wrk_nemo ! working arrays
14	USE diatmb
15	USE diurnal_bulk
16	USE cool_skin
17	USE ioipsl
18
19
20	IMPLICIT NONE
21	PRIVATE
22
23	LOGICAL , PUBLIC :: ln_diaopfoam !: Diaopfoam output
24	LOGICAL , PUBLIC :: ln_diaopfoam_Tzero !: Diaopfoam first time step output
25	PUBLIC dia_diaopfoam_init ! routine called by nemogcm.F90
26	PUBLIC dia_diaopfoam ! routine called by diawri.F90
27	PUBLIC calc_max_cur ! routine called by diaopfoam.F90
28
29	!! * Substitutions
30	# include "domzgr_substitute.h90"
31
32	!!----------------------------------------------------------------------
33	!! NEMO/OPA 3.6 , NEMO Consortium (2014)
34	!! $Id$
35	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
36	!!----------------------------------------------------------------------
37	CONTAINS
38
39	SUBROUTINE dia_diaopfoam_init
40	!!---------------------------------------------------------------------------
41	!! * ROUTINE dia_wri_diaop_init *
42	!!
43	!! ** Purpose: Initialization of diaopfoam namelist
44	!!
45	!! ** Method : Read namelist
46	!! History
47	!! 3.4 ! 03-14 (P. Sykes) Routine to initialize dia_wri_diaop
48	!!---------------------------------------------------------------------------
49	!!
50	INTEGER :: ios ! Local integer output status for namelist read
51	INTEGER :: ierror ! local integer
52	!!
53	NAMELIST/nam_diadiaop/ ln_diaopfoam,ln_diaopfoam_Tzero
54	!!
55	!!----------------------------------------------------------------------
56	!
57	ln_diaopfoam = .false. ! default value for diaopfoam stream
58	ln_diaopfoam_Tzero = .false. ! default value for diaopfoam Tzero stream
59	REWIND ( numnam_ref ) ! Read Namelist nam_diadiaop in reference namelist : 3D hourly diagnostics
60	READ ( numnam_ref, nam_diadiaop, IOSTAT=ios, ERR= 901 )
61	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in reference namelist', lwp )
62
63	REWIND( numnam_cfg ) ! Namelist nam_diadiaop in configuration namelist 3D hourly diagnostics
64	READ ( numnam_cfg, nam_diadiaop, IOSTAT = ios, ERR = 902 )
65	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in configuration namelist', lwp )
66	IF(lwm) WRITE ( numond, nam_diadiaop )
67	!
68	IF(lwp) THEN ! Control print
69	WRITE(numout,*)
70	WRITE(numout,*) 'dia_diaopfoam_init : Output Diaopfoam Diagnostics'
71	WRITE(numout,*) '~~~~~~~~~~~~'
72	WRITE(numout,*) ' Namelist nam_diadiaop : set diaopfoam outputs '
73	WRITE(numout,*) ' Switch for diaopfoam diagnostics (T) or not (F) ln_diaopfoam = ', ln_diaopfoam
74	WRITE(numout,*) ' Switch for diaopfoam first timestep diagnostics (T) or not (F) ln_diaopfoam_Tzero = ', ln_diaopfoam_Tzero
75	ENDIF
76	END SUBROUTINE dia_diaopfoam_init
77
78	SUBROUTINE dia_diaopfoam( kt )
79	!!----------------------------------------------------------------------
80	!! * ROUTINE dia_diaopfoam *
81	!! ** Purpose : Write 3D hourly diagnostics for operational use
82	!!
83	!!
84	!! History :
85	!! 3.6 ! 11-16 (P. Sykes)
86	!!
87	!!--------------------------------------------------------------------
88	IMPLICIT NONE
89
90	INTEGER, INTENT( in ) :: kt ! ocean time-step index
91
92	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
93	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
94	REAL(wp) :: zmdi
95	REAL(wp), POINTER, DIMENSION(:,:) :: zwu
96	REAL(wp), POINTER, DIMENSION(:,:) :: zwv
97	REAL(wp), POINTER, DIMENSION(:,:) :: zwz
98
99	CALL wrk_alloc( jpi , jpj , zwu )
100	CALL wrk_alloc( jpi , jpj , zwv )
101	CALL wrk_alloc( jpi , jpj , zwz )
102
103	zmdi=1.e+20 ! missing data indicator for masking
104
105	! Diaopfoam stream if needed
106	IF (ln_diaopfoam) THEN
107	IF ( kt .eq. nn_it000 .AND. ln_diaopfoam_Tzero ) THEN
108	IF(lwp) WRITE(numout,*) 'diaopfoam: writing T0 at kt = ', kt
109	CALL dia_diaopfoam_zero( 'Tzero', kt )
110	ENDIF
111
112	CALL theta2t
113	zw3d(:,:,:)=rinsitu_t(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
114	CALL iom_put( "insitut_op" , zw3d(:,:,:) ) ! insitu temperature
115	zw3d(:,:,:)=tsn(:,:,:,jp_tem)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
116	CALL iom_put( "toce_op" , zw3d(:,:,:) ) ! temperature
117	zw3d(:,:,:)=tsn(:,:,:,jp_sal)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
118	CALL iom_put( "soce_op" , zw3d(:,:,:) ) ! salinity
119	IF (ln_diurnal) THEN
120	CALL iom_put( "sst_wl_op" , x_dsst ) ! warm layer
121	CALL iom_put( "sst_cs_op" , x_csdsst ) ! cool skin
122	ENDIF
123	zw2d(:,:)=sshn(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
124	CALL iom_put( "ssh_op" , zw2d(:,:) ) ! sea surface height
125	zw3d(:,:,:)=un(:,:,:)umask(:,:,:) + zmdi(1.0-umask(:,:,:))
126	CALL iom_put( "uoce_op" , zw3d ) ! i-current
127	zw3d(:,:,:)=vn(:,:,:)vmask(:,:,:) + zmdi(1.0-vmask(:,:,:))
128	CALL iom_put( "voce_op" , zw3d ) ! j-current
129	!CALL iom_put( "woce_op" , wn ) ! k-current
130	CALL calc_max_cur(zwu,zwv,zwz,zmdi)
131	zw2d(:,:)=zwu(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
132	CALL iom_put( "maxu" , zw2d ) ! max u current
133	zw2d(:,:)=zwv(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
134	CALL iom_put( "maxv" , zw2d ) ! max v current
135	zw2d(:,:)=zwz(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
136	CALL iom_put( "maxz" , zw2d ) ! max current depth
137	ENDIF
138	END SUBROUTINE dia_diaopfoam
139
140	SUBROUTINE calc_max_cur(zmax_u, zmax_v, zmax_z, inmdi)
141	!!---------------------------------------------------------------------
142	!! * ROUTINE calc_max_cur *
143	!!
144	!! ** Purpose : To locate within the water column the magnitude and
145	!! vertical location of the strongest horizontal
146	!! current. The vertical component is ignored since it
147	!! is an order of magnitude smaller than the horizontal
148	!! flow, in general.
149	!!
150	!! ** Method : A. Map U,V to T-grid.
151	!! B. Calculate the magnitude of the current for every
152	!! grid cell.
153	!! C. Locate the vertical index using FORTRAN's builtin
154	!! MAXLOC function.
155	!! D. Copy the U,V,Z components of the relevant
156	!! indices.
157	!!
158	!! ** Returns : Value of u, v component and depth of maximum
159	!! horizontal current on T-grid.
160	!!
161	!!---------------------------------------------------------------------
162	IMPLICIT NONE
163	REAL(wp), DIMENSION(jpi, jpj), INTENT( OUT) :: zmax_u, zmax_v, zmax_z
164	REAL(wp), DIMENSION(jpi, jpj, jpk) :: zmax_u_t, zmax_v_t
165	REAL(wp), DIMENSION(jpi, jpj, jpk) :: zcmag
166	REAL(wp), DIMENSION(jpi, jpj) :: zmaxk
167	REAL(wp), INTENT(IN ) :: inmdi
168	INTEGER :: ji, jj, jk
169
170	! Initialise output arrays
171	zmax_u(:, :) = inmdi
172	zmax_v(:, :) = inmdi
173	zmax_z(:, :) = inmdi
174	! Map to T-grid
175	zmax_u_t(:, :, :) = 0._wp
176	zmax_v_t(:, :, :) = 0._wp
177	DO jk = 1,jpk
178	DO jj = 2,jpj
179	DO ji = 2,jpi
180	zmax_u_t(ji, jj, jk) = 0.5 * (un(ji, jj, jk) + un(ji-1, jj, jk))
181	zmax_v_t(ji, jj, jk) = 0.5 * (vn(ji, jj, jk) + vn(ji, jj-1, jk))
182	END DO
183	END DO
184	END DO
185	! Calculate absolute velocity
186	zcmag = sqrt((zmax_u_t)2 + (zmax_v_t)2)
187	! Find max. current
188	zmaxk = maxloc(zcmag, dim=3)
189	! Output values
190	DO jj = 1,jpj
191	DO ji = 1,jpi
192	zmax_u(ji, jj) = zmax_u_t(ji, jj, INT(zmaxk(ji, jj)))
193	zmax_v(ji, jj) = zmax_v_t(ji, jj, INT(zmaxk(ji, jj)))
194	zmax_z(ji, jj) = fsdept(ji, jj, INT(zmaxk(ji, jj)))
195	END DO
196	END DO
197	END SUBROUTINE calc_max_cur
198
199	SUBROUTINE dia_diaopfoam_zero( cdfile_name, kt )
200	!!---------------------------------------------------------------------
201	!! * ROUTINE dia_diaopfoam_zero *
202	!!
203	!! ** Purpose : create a NetCDF file named cdfile_name which contains
204	!! the instantaneous ocean state at the first tiome step.
205	!!
206	!! ** Method : NetCDF files using ioipsl
207	!!----------------------------------------------------------------------
208	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
209	INTEGER , INTENT( in ) :: kt ! ocean time-step index
210	!!
211	CHARACTER (len=32) :: clhstnam
212	CHARACTER (len=40) :: clop
213	INTEGER :: iimi, iima, ipk, ijmi, ijma ! local integers
214	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
215	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
216	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
217	INTEGER :: id_i , nz_i, nh_i
218	INTEGER, DIMENSION(1) :: idex ! local workspace
219	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
220	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
221	INTEGER :: ierr
222	INTEGER :: jkbot, jj, ji
223	REAL(wp) :: zsto, zout, zmax, zjulian, zdt
224	REAL(wp) :: zmdi
225	REAL(wp), POINTER, DIMENSION(:,:) :: zwu
226	REAL(wp), POINTER, DIMENSION(:,:) :: zwv
227	REAL(wp), POINTER, DIMENSION(:,:) :: zwz
228	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
229	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
230	REAL(wp), DIMENSION(jpi,jpj) :: z2d
231
232	!!----------------------------------------------------------------------
233
234	! -----------------------------------------------------------------
235	! 0. Allocations
236	! -----------------------------------------------------------------
237	ierr = 0
238	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
239	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
240	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr )
241	IF( lk_mpp ) CALL mpp_sum( ierr )
242	IF( ierr /= 0 ) THEN
243	CALL ctl_stop('dia_diaopfoam_zero: failed to allocate arrays')
244	RETURN
245	ENDIF
246	CALL wrk_alloc( jpi , jpj , zwu )
247	CALL wrk_alloc( jpi , jpj , zwv )
248	CALL wrk_alloc( jpi , jpj , zwz )
249
250	zmdi=1.e+20 ! missing data indicator for masking
251
252	! -----------------------------------------------------------------
253	! 1. Define NETCDF files and fields at beginning of first time step
254	! -----------------------------------------------------------------
255
256	! Define indices of the horizontal output zoom and vertical limit storage
257	iimi = 1 ; iima = jpi
258	ijmi = 1 ; ijma = jpj
259	ipk = jpk
260
261	! Define frequency of output and means
262	zdt = rdt
263	zsto = rdt
264	clop = "inst(x)" ! no use of the mask value (require less cpu time)
265	zout = rdt
266	zmax = ( nitend - nit000 + 1 ) * zdt
267
268	! Compute julian date from starting date of the run
269	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) ! time axis
270	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
271
272	! Define the T grid FILE ( nid_T )
273	clhstnam = TRIM(cdfile_name)//".grid_T"
274	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
275	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
276	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
277	& nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
278	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
279	& "m", ipk, gdept_1d, nz_T, "down" )
280	! ! Index of ocean points
281	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
282	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
283
284	! Define the U grid FILE ( nid_U )
285	clhstnam = TRIM(cdfile_name)//".grid_U"
286	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
287	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
288	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
289	& nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
290	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdepu
291	& "m", ipk, gdept_1d, nz_U, "down" )
292	! ! Index of ocean points
293	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
294	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
295
296	! Define the V grid FILE ( nid_V )
297	clhstnam = TRIM(cdfile_name)//".grid_V"
298	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
299	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
300	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
301	& nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
302	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdepv
303	& "m", ipk, gdept_1d, nz_V, "down" )
304	! ! Index of ocean points
305	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
306	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
307
308
309	! -----------------------------------------------------------------
310	! 2. Declare all the output fields as NETCDF variables
311	! -----------------------------------------------------------------
312
313	! !!! nid_T : 3D
314	!CALL histdef( nid_T, "votempis", "Insitu Temperature" , "C" , & !
315	! & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
316	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
317	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
318	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
319	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
320	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! sshn
321	& jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout )
322	CALL histdef( nid_T, "votempis", "Insitu Temperature" , "C" , & ! rinsitu_t
323	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
324	CALL histdef( nid_T, "maxu" , "Max Zonal Current" , "m/s" , & ! zwu
325	& jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout )
326	CALL histdef( nid_T, "maxv" , "Max Meridional Current" , "m/s" , & ! zwv
327	& jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout )
328	CALL histdef( nid_T, "maxz" , "Max Current Depth" , "m/s" , & ! zwz
329	& jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout )
330	CALL histdef( nid_T, "sbt" , "Bottom Temperature" , "C" , & ! sbt
331	& jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout )
332	CALL histend( nid_T, snc4chunks=snc4set )
333
334	! !!! nid_U : 3D
335	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
336	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
337	CALL histend( nid_U, snc4chunks=snc4set )
338
339	! !!! nid_V : 3D
340	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
341	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
342	CALL histend( nid_V, snc4chunks=snc4set )
343
344
345	! -----------------------------------------------------------------
346	! 3. Write the data
347	! -----------------------------------------------------------------
348
349	idex(1) = 1
350	zw3d(:,:,:)=tsn(:,:,:,jp_tem)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
351	CALL histwrite( nid_T, "votemper", kt, zw3d(:,:,:) , jpijpjjpk, idex ) ! now temperature
352	zw3d(:,:,:)=tsn(:,:,:,jp_sal)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
353	CALL histwrite( nid_T, "vosaline", kt, zw3d(:,:,:) , jpijpjjpk, idex ) ! now salinity
354	zw2d(:,:)=sshn(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
355	CALL histwrite( nid_T, "sossheig", kt, zw2d(:,:) , jpi*jpj , idex ) ! sea surface height
356	CALL theta2t
357	zw3d(:,:,:)=rinsitu_t(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
358	CALL histwrite( nid_T, "votempis", kt, zw3d(:,:,:) , jpijpjjpk, idex ) ! now insitu-temperature
359	CALL calc_max_cur(zwu,zwv,zwz,zmdi)
360	CALL lbc_lnk( zwu, 'T', 1. )
361	CALL lbc_lnk( zwv, 'T', 1. )
362	CALL lbc_lnk( zwz, 'T', 1. )
363	zw2d(:,:)=zwu(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
364	CALL histwrite( nid_T, "maxu" , kt, zw2d , jpi*jpj , idex ) ! max u-current
365	zw2d(:,:)=zwv(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
366	CALL histwrite( nid_T, "maxv" , kt, zw2d , jpi*jpj , idex ) ! max v-current
367	zw2d(:,:)=zwz(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
368	CALL histwrite( nid_T, "maxz" , kt, zw2d , jpi*jpj , idex ) ! max current depth
369	DO jj = 1, jpj
370	DO ji = 1, jpi
371	jkbot = mbkt(ji,jj)
372	z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)tmask(ji,jj,jkbot) + zmdi(1.0-tmask(ji,jj,jkbot))
373	END DO
374	END DO
375	CALL histwrite( nid_T, "sbt" , kt, z2d , jpi*jpj , idex ) ! sbt
376	! U file
377	zw3d(:,:,:)=un(:,:,:)umask(:,:,:) + zmdi(1.0-umask(:,:,:))
378	CALL histwrite( nid_U, "vozocrtx", kt, zw3d(:,:,:) , jpijpjjpk, idex ) ! now i-velocity
379	! V file
380	zw3d(:,:,:)=vn(:,:,:)vmask(:,:,:) + zmdi(1.0-vmask(:,:,:))
381	CALL histwrite( nid_V, "vomecrty", kt, zw3d(:,:,:) , jpijpjjpk, idex ) ! now j-velocity
382
383
384	! -----------------------------------------------------------------
385	! 4. Close the files
386	! -----------------------------------------------------------------
387
388	CALL histclo( nid_T )
389	CALL histclo( nid_U )
390	CALL histclo( nid_V )
391
392	END SUBROUTINE dia_diaopfoam_zero
393
394	END MODULE diaopfoam

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_PS44/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90

Download in other formats: