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.

dia25h.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
File size: 17.2 KB

Line
1	MODULE dia25h
2	!!======================================================================
3	!! * MODULE diaharm *
4	!! Harmonic analysis of tidal constituents
5	!!======================================================================
6	!! History : 3.6 ! 2014 (E O'Dea) Original code
7	!!----------------------------------------------------------------------
8	USE oce ! ocean dynamics and tracers variables
9	USE dom_oce ! ocean space and time domain
10	USE in_out_manager ! I/O units
11	USE iom ! I/0 library
12	USE wrk_nemo ! working arrays
13	#if defined key_zdftke
14	USE zdf_oce, ONLY: en
15	#endif
16	USE zdf_oce, ONLY: avt, avm
17	#if defined key_zdfgls
18	USE zdf_oce, ONLY: en
19	USE zdfgls, ONLY: mxln
20	#endif
21
22	IMPLICIT NONE
23	PRIVATE
24
25	LOGICAL , PUBLIC :: ln_dia25h !: 25h mean output
26	PUBLIC dia_25h_init ! routine called by nemogcm.F90
27	PUBLIC dia_25h ! routine called by diawri.F90
28
29	!! * variables for calculating 25-hourly means
30	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tn_25h , sn_25h
31	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:) :: sshn_25h
32	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: un_25h , vn_25h , wn_25h
33	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avt_25h , avm_25h
34	#if defined key_zdfgls \|\| key_zdftke
35	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: en_25h
36	#endif
37	#if defined key_zdfgls
38	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: rmxln_25h
39	#endif
40	INTEGER, SAVE :: cnt_25h ! Counter for 25 hour means
41
42
43
44	!!----------------------------------------------------------------------
45	!! NEMO/OPA 3.6 , NEMO Consortium (2014)
46	!! $Id:$
47	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
48	!!----------------------------------------------------------------------
49	CONTAINS
50
51	SUBROUTINE dia_25h_init
52	!!---------------------------------------------------------------------------
53	!! * ROUTINE dia_25h_init *
54	!!
55	!! ** Purpose: Initialization of 25h mean namelist
56	!!
57	!! ** Method : Read namelist
58	!! History
59	!! 3.6 ! 08-14 (E. O'Dea) Routine to initialize dia_25h
60	!!---------------------------------------------------------------------------
61	!!
62	INTEGER :: ios ! Local integer output status for namelist read
63	INTEGER :: ierror ! Local integer for memory allocation
64	INTEGER :: ji, jj, jk
65	!
66	NAMELIST/nam_dia25h/ ln_dia25h
67	!!----------------------------------------------------------------------
68	!
69	REWIND ( numnam_ref ) ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
70	READ ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
71	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp )
72
73	REWIND( numnam_cfg ) ! Namelist nam_dia25h in configuration namelist 25hour diagnostics
74	READ ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
75	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp )
76	IF(lwm) WRITE ( numond, nam_dia25h )
77
78	IF(lwp) THEN ! Control print
79	WRITE(numout,*)
80	WRITE(numout,*) 'dia_25h_init : Output 25 hour mean diagnostics'
81	WRITE(numout,*) '~~~~~~~~~~~~'
82	WRITE(numout,*) 'Namelist nam_dia25h : set 25h outputs '
83	WRITE(numout,*) 'Switch for 25h diagnostics (T) or not (F) ln_dia25h = ', ln_dia25h
84	ENDIF
85	IF( .NOT. ln_dia25h ) RETURN
86	! ------------------- !
87	! 1 - Allocate memory !
88	! ------------------- !
89	ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror )
90	IF( ierror > 0 ) THEN
91	CALL ctl_stop( 'dia_25h: unable to allocate tn_25h' ) ; RETURN
92	ENDIF
93	ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror )
94	IF( ierror > 0 ) THEN
95	CALL ctl_stop( 'dia_25h: unable to allocate sn_25h' ) ; RETURN
96	ENDIF
97	ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror )
98	IF( ierror > 0 ) THEN
99	CALL ctl_stop( 'dia_25h: unable to allocate un_25h' ) ; RETURN
100	ENDIF
101	ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror )
102	IF( ierror > 0 ) THEN
103	CALL ctl_stop( 'dia_25h: unable to allocate vn_25h' ) ; RETURN
104	ENDIF
105	ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror )
106	IF( ierror > 0 ) THEN
107	CALL ctl_stop( 'dia_25h: unable to allocate wn_25h' ) ; RETURN
108	ENDIF
109	ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror )
110	IF( ierror > 0 ) THEN
111	CALL ctl_stop( 'dia_25h: unable to allocate avt_25h' ) ; RETURN
112	ENDIF
113	ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror )
114	IF( ierror > 0 ) THEN
115	CALL ctl_stop( 'dia_25h: unable to allocate avm_25h' ) ; RETURN
116	ENDIF
117	# if defined key_zdfgls \|\| defined key_zdftke
118	ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
119	IF( ierror > 0 ) THEN
120	CALL ctl_stop( 'dia_25h: unable to allocate en_25h' ) ; RETURN
121	ENDIF
122	#endif
123	# if defined key_zdfgls
124	ALLOCATE( rmxln_25h(jpi,jpj,jpk), STAT=ierror )
125	IF( ierror > 0 ) THEN
126	CALL ctl_stop( 'dia_25h: unable to allocate rmxln_25h' ) ; RETURN
127	ENDIF
128	#endif
129	ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror )
130	IF( ierror > 0 ) THEN
131	CALL ctl_stop( 'dia_25h: unable to allocate sshn_25h' ) ; RETURN
132	ENDIF
133	! ------------------------- !
134	! 2 - Assign Initial Values !
135	! ------------------------- !
136	cnt_25h = 1 ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)
137	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
138	DO jk = 1, jpk
139	DO jj = 1, jpj
140	DO ji = 1, jpi
141	tn_25h(ji,jj,jk) = tsb(ji,jj,jk,jp_tem)
142	sn_25h(ji,jj,jk) = tsb(ji,jj,jk,jp_sal)
143	sshn_25h(ji,jj) = sshb(ji,jj)
144	un_25h(ji,jj,jk) = ub(ji,jj,jk)
145	vn_25h(ji,jj,jk) = vb(ji,jj,jk)
146	wn_25h(ji,jj,jk) = wn(ji,jj,jk)
147	avt_25h(ji,jj,jk) = avt(ji,jj,jk)
148	avm_25h(ji,jj,jk) = avm(ji,jj,jk)
149	# if defined key_zdfgls \|\| defined key_zdftke
150	en_25h(ji,jj,jk) = en(ji,jj,jk)
151	#endif
152	# if defined key_zdfgls
153	rmxln_25h(ji,jj,jk) = mxln(ji,jj,jk)
154	#endif
155	END DO
156	END DO
157	END DO
158	#if defined key_lim3 \|\| defined key_lim2
159	CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
160	#endif
161
162	! -------------------------- !
163	! 3 - Return to dia_wri !
164	! -------------------------- !
165
166
167	END SUBROUTINE dia_25h_init
168
169
170	SUBROUTINE dia_25h( kt )
171	!!----------------------------------------------------------------------
172	!! * ROUTINE dia_25h *
173	!!
174	!!
175	!!--------------------------------------------------------------------
176	!!
177	!! ** Purpose : Write diagnostics with M2/S2 tide removed
178	!!
179	!! ** Method :
180	!! 25hr mean outputs for shelf seas
181	!!
182	!! History :
183	!! ?.0 ! 07-04 (A. Hines) New routine, developed from dia_wri_foam
184	!! 3.4 ! 02-13 (J. Siddorn) Routine taken from old dia_wri_foam
185	!! 3.6 ! 08-14 (E. O'Dea) adapted for VN3.6
186	!!----------------------------------------------------------------------
187	!! * Modules used
188
189	IMPLICIT NONE
190
191	!! * Arguments
192	INTEGER, INTENT( in ) :: kt ! ocean time-step index
193
194
195	!! * Local declarations
196	INTEGER :: ji, jj, jk
197
198	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
199	REAL(wp) :: zsto, zout, zmax, zjulian, zmdi ! temporary reals
200	INTEGER :: i_steps ! no of timesteps per hour
201	REAL(wp), DIMENSION(jpi,jpj ) :: zw2d, un_dm, vn_dm ! temporary workspace
202	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! temporary workspace
203	REAL(wp), DIMENSION(jpi,jpj,3) :: zwtmb ! temporary workspace
204	INTEGER :: iyear0, nimonth0,iday0 ! start year,imonth,day
205
206	!!----------------------------------------------------------------------
207
208	! 0. Initialisation
209	! -----------------
210	! Define frequency of summing to create 25 h mean
211	IF( MOD( 3600,INT(rdt) ) == 0 ) THEN
212	i_steps = 3600/INT(rdt)
213	ELSE
214	CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
215	ENDIF
216
217	#if defined key_lim3 \|\| defined key_lim2
218	CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice')
219	#endif
220
221	! local variable for debugging
222	ll_print = ll_print .AND. lwp
223
224	! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours
225	! every day
226	IF( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 ) THEN
227
228	IF (lwp) THEN
229	WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt
230	WRITE(numout,*) '~~~~~~~~~~~~ '
231	ENDIF
232
233	!$OMP PARALLEL
234	!$OMP DO schedule(static) private(jj, ji)
235	DO jj = 1, jpj
236	DO ji = 1, jpi
237	sshn_25h(ji,jj) = sshn_25h(ji,jj) + sshn (ji,jj)
238	END DO
239	END DO
240	!$OMP END DO NOWAIT
241	!$OMP DO schedule(static) private(jk, jj, ji)
242	DO jk = 1, jpk
243	DO jj = 1, jpj
244	DO ji = 1, jpi
245	tn_25h(ji,jj,jk) = tn_25h(ji,jj,jk) + tsn(ji,jj,jk,jp_tem)
246	sn_25h(ji,jj,jk) = sn_25h(ji,jj,jk) + tsn(ji,jj,jk,jp_sal)
247	un_25h(ji,jj,jk) = un_25h(ji,jj,jk) + un(ji,jj,jk)
248	vn_25h(ji,jj,jk) = vn_25h(ji,jj,jk) + vn(ji,jj,jk)
249	wn_25h(ji,jj,jk) = wn_25h(ji,jj,jk) + wn(ji,jj,jk)
250	avt_25h(ji,jj,jk) = avt_25h(ji,jj,jk) + avt(ji,jj,jk)
251	avm_25h(ji,jj,jk) = avm_25h(ji,jj,jk) + avm(ji,jj,jk)
252	# if defined key_zdfgls \|\| defined key_zdftke
253	en_25h(ji,jj,jk) = en_25h(ji,jj,jk) + en(ji,jj,jk)
254	#endif
255	# if defined key_zdfgls
256	rmxln_25h(ji,jj,jk) = rmxln_25h(ji,jj,jk) + mxln(ji,jj,jk)
257	#endif
258	END DO
259	END DO
260	END DO
261	!$OMP END PARALLEL
262	cnt_25h = cnt_25h + 1
263
264	IF (lwp) THEN
265	WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
266	ENDIF
267
268	ENDIF ! MOD( kt, i_steps ) == 0
269
270	! Write data for 25 hour mean output streams
271	IF( cnt_25h .EQ. 25 .AND. MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
272
273	IF(lwp) THEN
274	WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
275	WRITE(numout,*) '~~~~~~~~~~~~ '
276	ENDIF
277
278	!$OMP PARALLEL
279	!$OMP DO schedule(static) private(jj, ji)
280	DO jj = 1, jpj
281	DO ji = 1, jpi
282	sshn_25h(ji,jj) = sshn_25h(ji,jj) / 25.0_wp
283	END DO
284	END DO
285	!$OMP END DO NOWAIT
286	!$OMP DO schedule(static) private(jk, jj, ji)
287	DO jk = 1, jpk
288	DO jj = 1, jpj
289	DO ji = 1, jpi
290	tn_25h(ji,jj,jk) = tn_25h(ji,jj,jk) / 25.0_wp
291	sn_25h(ji,jj,jk) = sn_25h(ji,jj,jk) / 25.0_wp
292	un_25h(ji,jj,jk) = un_25h(ji,jj,jk) / 25.0_wp
293	vn_25h(ji,jj,jk) = vn_25h(ji,jj,jk) / 25.0_wp
294	wn_25h(ji,jj,jk) = wn_25h(ji,jj,jk) / 25.0_wp
295	avt_25h(ji,jj,jk) = avt_25h(ji,jj,jk) / 25.0_wp
296	avm_25h(ji,jj,jk) = avm_25h(ji,jj,jk) / 25.0_wp
297	# if defined key_zdfgls \|\| defined key_zdftke
298	en_25h(ji,jj,jk) = en_25h(ji,jj,jk) / 25.0_wp
299	#endif
300	# if defined key_zdfgls
301	rmxln_25h(ji,jj,jk) = rmxln_25h(ji,jj,jk) / 25.0_wp
302	#endif
303	END DO
304	END DO
305	END DO
306	!$OMP END PARALLEL
307
308	IF (lwp) WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
309	zmdi=1.e+20 !missing data indicator for masking
310	! write tracers (instantaneous)
311	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
312	DO jk = 1, jpk
313	DO jj = 1, jpj
314	DO ji = 1, jpi
315	zw3d(ji,jj,jk) = tn_25h(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
316	END DO
317	END DO
318	END DO
319	CALL iom_put("temper25h", zw3d) ! potential temperature
320	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
321	DO jk = 1, jpk
322	DO jj = 1, jpj
323	DO ji = 1, jpi
324	zw3d(ji,jj,jk) = sn_25h(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
325	END DO
326	END DO
327	END DO
328	CALL iom_put( "salin25h", zw3d ) ! salinity
329	!$OMP PARALLEL DO schedule(static) private(jj, ji)
330	DO jj = 1, jpj
331	DO ji = 1, jpi
332	zw2d(ji,jj) = sshn_25h(ji,jj)tmask(ji,jj,1) + zmdi(1.0-tmask(ji,jj,1))
333	END DO
334	END DO
335	CALL iom_put( "ssh25h", zw2d ) ! sea surface
336
337
338	! Write velocities (instantaneous)
339	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
340	DO jk = 1, jpk
341	DO jj = 1, jpj
342	DO ji = 1, jpi
343	zw3d(ji,jj,jk) = un_25h(ji,jj,jk)umask(ji,jj,jk) + zmdi(1.0-umask(ji,jj,jk))
344	END DO
345	END DO
346	END DO
347	CALL iom_put("vozocrtx25h", zw3d) ! i-current
348	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
349	DO jk = 1, jpk
350	DO jj = 1, jpj
351	DO ji = 1, jpi
352	zw3d(ji,jj,jk) = vn_25h(ji,jj,jk)vmask(ji,jj,jk) + zmdi(1.0-vmask(ji,jj,jk))
353	END DO
354	END DO
355	END DO
356	CALL iom_put("vomecrty25h", zw3d ) ! j-current
357	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
358	DO jk = 1, jpk
359	DO jj = 1, jpj
360	DO ji = 1, jpi
361	zw3d(ji,jj,jk) = wn_25h(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
362	END DO
363	END DO
364	END DO
365	CALL iom_put("vomecrtz25h", zw3d ) ! k-current
366	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
367	DO jk = 1, jpk
368	DO jj = 1, jpj
369	DO ji = 1, jpi
370	zw3d(ji,jj,jk) = avt_25h(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
371	END DO
372	END DO
373	END DO
374	CALL iom_put("avt25h", zw3d ) ! diffusivity
375	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
376	DO jk = 1, jpk
377	DO jj = 1, jpj
378	DO ji = 1, jpi
379	zw3d(ji,jj,jk) = avm_25h(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
380	END DO
381	END DO
382	END DO
383	CALL iom_put("avm25h", zw3d) ! viscosity
384	#if defined key_zdftke \|\| defined key_zdfgls
385	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
386	DO jk = 1, jpk
387	DO jj = 1, jpj
388	DO ji = 1, jpi
389	zw3d(ji,jj,jk) = en_25h(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
390	END DO
391	END DO
392	END DO
393	CALL iom_put("tke25h", zw3d) ! tke
394	#endif
395	#if defined key_zdfgls
396	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
397	DO jk = 1, jpk
398	DO jj = 1, jpj
399	DO ji = 1, jpi
400	zw3d(ji,jj,jk) = rmxln_25h(ji,jj,jk)tmask(ji,jj,jk) + zmdi(1.0-tmask(ji,jj,jk))
401	END DO
402	END DO
403	END DO
404	CALL iom_put( "mxln25h",zw3d)
405	#endif
406
407	! After the write reset the values to cnt=1 and sum values equal current value
408	!$OMP PARALLEL
409	!$OMP DO schedule(static) private(jj, ji)
410	DO jj = 1, jpj
411	DO ji = 1, jpi
412	sshn_25h(ji,jj) = sshn (ji,jj)
413	END DO
414	END DO
415	!$OMP END DO NOWAIT
416	!$OMP DO schedule(static) private(jk, jj, ji)
417	DO jk = 1, jpk
418	DO jj = 1, jpj
419	DO ji = 1, jpi
420	tn_25h(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
421	sn_25h(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
422	un_25h(ji,jj,jk) = un(ji,jj,jk)
423	vn_25h(ji,jj,jk) = vn(ji,jj,jk)
424	wn_25h(ji,jj,jk) = wn(ji,jj,jk)
425	avt_25h(ji,jj,jk) = avt(ji,jj,jk)
426	avm_25h(ji,jj,jk) = avm(ji,jj,jk)
427	# if defined key_zdfgls \|\| defined key_zdftke
428	en_25h(ji,jj,jk) = en(ji,jj,jk)
429	#endif
430	# if defined key_zdfgls
431	rmxln_25h(ji,jj,jk) = mxln(ji,jj,jk)
432	#endif
433	END DO
434	END DO
435	END DO
436	!$OMP END PARALLEL
437	cnt_25h = 1
438	IF (lwp) WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
439
440	ENDIF ! cnt_25h .EQ. 25 .AND. MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
441
442
443	END SUBROUTINE dia_25h
444
445	!!======================================================================
446	END MODULE dia25h

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

Download in other formats: