Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

lib_fortran.F90

Last change on this file was 15376, checked in by clem, 3 years ago
add functions glob_min_vec and glob_max_vec. They would need to be merged into the generic module lib_fortran_generic.h90 at some point, along with glob_sum_vec
Property svn:keywords set to `Id`
File size: 32.6 KB

Line
1	MODULE lib_fortran
2	!!======================================================================
3	!! * MODULE lib_fortran *
4	!! Fortran utilities: includes some low levels fortran functionality
5	!!======================================================================
6	!! History : 3.2 ! 2010-05 (M. Dunphy, R. Benshila) Original code
7	!! 3.4 ! 2013-06 (C. Rousset) add glob_min, glob_max
8	!! + 3d dim. of input is fexible (jpk, jpl...)
9	!! 4.0 ! 2016-06 (T. Lovato) double precision global sum by default
10	!!----------------------------------------------------------------------
11
12	!!----------------------------------------------------------------------
13	!! glob_sum : generic interface for global masked summation over
14	!! the interior domain for 1 or 2 2D or 3D arrays
15	!! it works only for T points
16	!! SIGN : generic interface for SIGN to overwrite f95 behaviour
17	!! of intrinsinc sign function
18	!!----------------------------------------------------------------------
19	USE par_oce ! Ocean parameter
20	USE dom_oce ! ocean domain
21	USE in_out_manager ! I/O manager
22	USE lib_mpp ! distributed memory computing
23	USE lbclnk ! ocean lateral boundary conditions
24
25	IMPLICIT NONE
26	PRIVATE
27
28	PUBLIC glob_sum ! used in many places (masked with tmask_i = ssmask * tmask_h)
29	PUBLIC glob_sum_full ! used in many places (masked with tmask_h, excluding all duplicated points halos+periodicity)
30	PUBLIC local_sum ! used in trcrad, local operation before glob_sum_delay
31	PUBLIC sum3x3 ! used in trcrad, do a sum over 3x3 boxes
32	PUBLIC DDPDD ! also used in closea module
33	PUBLIC glob_min, glob_max
34	PUBLIC glob_sum_vec
35	PUBLIC glob_sum_full_vec
36	PUBLIC glob_min_vec, glob_max_vec
37	#if defined key_nosignedzero
38	PUBLIC SIGN
39	#endif
40
41	INTERFACE glob_sum
42	MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d
43	END INTERFACE
44	INTERFACE glob_sum_full
45	MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d
46	END INTERFACE
47	INTERFACE local_sum
48	MODULE PROCEDURE local_sum_2d, local_sum_3d
49	END INTERFACE
50	INTERFACE sum3x3
51	MODULE PROCEDURE sum3x3_2d, sum3x3_3d
52	END INTERFACE
53	INTERFACE glob_min
54	MODULE PROCEDURE glob_min_2d, glob_min_3d
55	END INTERFACE
56	INTERFACE glob_max
57	MODULE PROCEDURE glob_max_2d, glob_max_3d
58	END INTERFACE
59	INTERFACE glob_sum_vec
60	MODULE PROCEDURE glob_sum_vec_3d, glob_sum_vec_4d
61	END INTERFACE
62	INTERFACE glob_sum_full_vec
63	MODULE PROCEDURE glob_sum_full_vec_3d, glob_sum_full_vec_4d
64	END INTERFACE
65	INTERFACE glob_min_vec
66	MODULE PROCEDURE glob_min_vec_3d, glob_min_vec_4d
67	END INTERFACE
68	INTERFACE glob_max_vec
69	MODULE PROCEDURE glob_max_vec_3d, glob_max_vec_4d
70	END INTERFACE
71
72	#if defined key_nosignedzero
73	INTERFACE SIGN
74	MODULE PROCEDURE SIGN_SCALAR, SIGN_ARRAY_1D, SIGN_ARRAY_2D, SIGN_ARRAY_3D, &
75	& SIGN_ARRAY_1D_A, SIGN_ARRAY_2D_A, SIGN_ARRAY_3D_A, &
76	& SIGN_ARRAY_1D_B, SIGN_ARRAY_2D_B, SIGN_ARRAY_3D_B
77	END INTERFACE
78	#endif
79
80	!! * Substitutions
81	# include "do_loop_substitute.h90"
82	!!----------------------------------------------------------------------
83	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
84	!! $Id$
85	!! Software governed by the CeCILL license (see ./LICENSE)
86	!!----------------------------------------------------------------------
87	CONTAINS
88
89	# define GLOBSUM_CODE
90
91	# define DIM_1d
92	# define FUNCTION_GLOBSUM glob_sum_1d
93	# include "lib_fortran_generic.h90"
94	# undef FUNCTION_GLOBSUM
95	# undef DIM_1d
96
97	# define DIM_2d
98	# define OPERATION_GLOBSUM
99	# define FUNCTION_GLOBSUM glob_sum_2d
100	# include "lib_fortran_generic.h90"
101	# undef FUNCTION_GLOBSUM
102	# undef OPERATION_GLOBSUM
103	# define OPERATION_FULL_GLOBSUM
104	# define FUNCTION_GLOBSUM glob_sum_full_2d
105	# include "lib_fortran_generic.h90"
106	# undef FUNCTION_GLOBSUM
107	# undef OPERATION_FULL_GLOBSUM
108	# undef DIM_2d
109
110	# define DIM_3d
111	# define OPERATION_GLOBSUM
112	# define FUNCTION_GLOBSUM glob_sum_3d
113	# include "lib_fortran_generic.h90"
114	# undef FUNCTION_GLOBSUM
115	# undef OPERATION_GLOBSUM
116	# define OPERATION_FULL_GLOBSUM
117	# define FUNCTION_GLOBSUM glob_sum_full_3d
118	# include "lib_fortran_generic.h90"
119	# undef FUNCTION_GLOBSUM
120	# undef OPERATION_FULL_GLOBSUM
121	# undef DIM_3d
122
123	# undef GLOBSUM_CODE
124
125
126	# define GLOBMINMAX_CODE
127
128	# define DIM_2d
129	# define OPERATION_GLOBMIN
130	# define FUNCTION_GLOBMINMAX glob_min_2d
131	# include "lib_fortran_generic.h90"
132	# undef FUNCTION_GLOBMINMAX
133	# undef OPERATION_GLOBMIN
134	# define OPERATION_GLOBMAX
135	# define FUNCTION_GLOBMINMAX glob_max_2d
136	# include "lib_fortran_generic.h90"
137	# undef FUNCTION_GLOBMINMAX
138	# undef OPERATION_GLOBMAX
139	# undef DIM_2d
140
141	# define DIM_3d
142	# define OPERATION_GLOBMIN
143	# define FUNCTION_GLOBMINMAX glob_min_3d
144	# include "lib_fortran_generic.h90"
145	# undef FUNCTION_GLOBMINMAX
146	# undef OPERATION_GLOBMIN
147	# define OPERATION_GLOBMAX
148	# define FUNCTION_GLOBMINMAX glob_max_3d
149	# include "lib_fortran_generic.h90"
150	# undef FUNCTION_GLOBMINMAX
151	# undef OPERATION_GLOBMAX
152	# undef DIM_3d
153	# undef GLOBMINMAX_CODE
154
155	! ! FUNCTION local_sum !
156
157	FUNCTION local_sum_2d( ptab )
158	!!----------------------------------------------------------------------
159	REAL(wp), INTENT(in ) :: ptab(:,:) ! array on which operation is applied
160	COMPLEX(dp) :: local_sum_2d
161	!
162	!!-----------------------------------------------------------------------
163	!
164	COMPLEX(dp):: ctmp
165	REAL(wp) :: ztmp
166	INTEGER :: ji, jj ! dummy loop indices
167	INTEGER :: ipi, ipj ! dimensions
168	!!-----------------------------------------------------------------------
169	!
170	ipi = SIZE(ptab,1) ! 1st dimension
171	ipj = SIZE(ptab,2) ! 2nd dimension
172	!
173	ctmp = CMPLX( 0.e0, 0.e0, wp ) ! warning ctmp is cumulated
174
175	DO jj = 1, ipj
176	DO ji = 1, ipi
177	ztmp = ptab(ji,jj) * tmask_i(ji,jj)
178	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
179	END DO
180	END DO
181	!
182	local_sum_2d = ctmp
183
184	END FUNCTION local_sum_2d
185
186	FUNCTION local_sum_3d( ptab )
187	!!----------------------------------------------------------------------
188	REAL(wp), INTENT(in ) :: ptab(:,:,:) ! array on which operation is applied
189	COMPLEX(dp) :: local_sum_3d
190	!
191	!!-----------------------------------------------------------------------
192	!
193	COMPLEX(dp):: ctmp
194	REAL(wp) :: ztmp
195	INTEGER :: ji, jj, jk ! dummy loop indices
196	INTEGER :: ipi, ipj, ipk ! dimensions
197	!!-----------------------------------------------------------------------
198	!
199	ipi = SIZE(ptab,1) ! 1st dimension
200	ipj = SIZE(ptab,2) ! 2nd dimension
201	ipk = SIZE(ptab,3) ! 3rd dimension
202	!
203	ctmp = CMPLX( 0.e0, 0.e0, wp ) ! warning ctmp is cumulated
204
205	DO jk = 1, ipk
206	DO jj = 1, ipj
207	DO ji = 1, ipi
208	ztmp = ptab(ji,jj,jk) * tmask_i(ji,jj)
209	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
210	END DO
211	END DO
212	END DO
213	!
214	local_sum_3d = ctmp
215
216	END FUNCTION local_sum_3d
217
218	! ! FUNCTION sum3x3 !
219
220	SUBROUTINE sum3x3_2d( p2d )
221	!!-----------------------------------------------------------------------
222	!! * routine sum3x3_2d *
223	!!
224	!! ** Purpose : sum over 3x3 boxes
225	!!----------------------------------------------------------------------
226	REAL(wp), DIMENSION (:,:), INTENT(inout) :: p2d
227	!
228	INTEGER :: ji, ji2, jj, jj2 ! dummy loop indices
229	!!----------------------------------------------------------------------
230	!
231	IF( SIZE(p2d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the first dimension is not equal to jpi' )
232	IF( SIZE(p2d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the second dimension is not equal to jpj' )
233	!
234	! work over the whole domain (guarantees all internal cells are set when nn_hls=2)
235	!
236	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
237	IF( MOD(mig(ji), 3) == MOD(nn_hls, 3) .AND. & ! 1st bottom left corner always at (Nis0-1, Njs0-1)
238	& MOD(mjg(jj), 3) == MOD(nn_hls, 3) ) THEN ! bottom left corner of a 3x3 box
239	ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1 ! right position of the box
240	jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1 ! upper position of the box
241	IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN ! the box is fully included in the local mpi domain
242	p2d(ji:ji2,jj:jj2) = SUM(p2d(ji:ji2,jj:jj2))
243	ENDIF
244	ENDIF
245	END_2D
246	CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1.0_wp )
247	! no need for 2nd exchange when nn_hls > 1
248	IF( nn_hls == 1 ) THEN
249	IF( mpiRnei(nn_hls,jpwe) > -1 ) THEN ! 1st column was changed during the previous call to lbc_lnk
250	IF( MOD(mig( 1), 3) == 1 ) & ! 1st box start at i=1 -> column 1 to 3 correctly computed locally
251	p2d( 1,:) = p2d( 2,:) ! previous lbc_lnk corrupted column 1 -> put it back using column 2
252	IF( MOD(mig( 1), 3) == 2 ) & ! 1st box start at i=3 -> column 1 and 2 correctly computed on west neighbourh
253	p2d( 2,:) = p2d( 1,:) ! previous lbc_lnk fix column 1 -> copy it to column 2
254	ENDIF
255	IF( mpiRnei(nn_hls,jpea) > -1 ) THEN
256	IF( MOD(mig(jpi-2), 3) == 1 ) p2d( jpi,:) = p2d(jpi-1,:)
257	IF( MOD(mig(jpi-2), 3) == 0 ) p2d(jpi-1,:) = p2d( jpi,:)
258	ENDIF
259	IF( mpiRnei(nn_hls,jpso) > -1 ) THEN
260	IF( MOD(mjg( 1), 3) == 1 ) p2d(:, 1) = p2d(:, 2)
261	IF( MOD(mjg( 1), 3) == 2 ) p2d(:, 2) = p2d(:, 1)
262	ENDIF
263	IF( mpiRnei(nn_hls,jpno) > -1 ) THEN
264	IF( MOD(mjg(jpj-2), 3) == 1 ) p2d(:, jpj) = p2d(:,jpj-1)
265	IF( MOD(mjg(jpj-2), 3) == 0 ) p2d(:,jpj-1) = p2d(:, jpj)
266	ENDIF
267	CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1.0_wp )
268	ENDIF
269
270	END SUBROUTINE sum3x3_2d
271
272	SUBROUTINE sum3x3_3d( p3d )
273	!!-----------------------------------------------------------------------
274	!! * routine sum3x3_3d *
275	!!
276	!! ** Purpose : sum over 3x3 boxes
277	!!----------------------------------------------------------------------
278	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: p3d
279	!
280	INTEGER :: ji, ji2, jj, jj2, jn ! dummy loop indices
281	INTEGER :: ipn ! Third dimension size
282	!!----------------------------------------------------------------------
283	!
284	IF( SIZE(p3d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the first dimension is not equal to jpi' )
285	IF( SIZE(p3d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the second dimension is not equal to jpj' )
286	ipn = SIZE(p3d,3)
287	!
288	DO jn = 1, ipn
289	!
290	! work over the whole domain (guarantees all internal cells are set when nn_hls=2)
291	!
292	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
293	IF( MOD(mig(ji), 3) == MOD(nn_hls, 3) .AND. & ! 1st bottom left corner always at (Nis0-1, Njs0-1)
294	& MOD(mjg(jj), 3) == MOD(nn_hls, 3) ) THEN ! bottom left corner of a 3x3 box
295	ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1 ! right position of the box
296	jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1 ! upper position of the box
297	IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN ! the box is fully included in the local mpi domain
298	p3d(ji:ji2,jj:jj2,jn) = SUM(p3d(ji:ji2,jj:jj2,jn))
299	ENDIF
300	ENDIF
301	END_2D
302	END DO
303	CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1.0_wp )
304	! no need for 2nd exchange when nn_hls > 1
305	IF( nn_hls == 1 ) THEN
306	IF( mpiRnei(nn_hls,jpwe) > -1 ) THEN ! 1st column was changed during the previous call to lbc_lnk
307	IF( MOD(mig( 1), 3) == 1 ) & ! 1st box start at i=1 -> column 1 to 3 correctly computed locally
308	p3d( 1,:,:) = p3d( 2,:,:) ! previous lbc_lnk corrupted column 1 -> put it back using column 2
309	IF( MOD(mig( 1), 3) == 2 ) & ! 1st box start at i=3 -> column 1 and 2 correctly computed on west neighbourh
310	p3d( 2,:,:) = p3d( 1,:,:) ! previous lbc_lnk fix column 1 -> copy it to column 2
311	ENDIF
312	IF( mpiRnei(nn_hls,jpea) > -1 ) THEN
313	IF( MOD(mig(jpi-2), 3) == 1 ) p3d( jpi,:,:) = p3d(jpi-1,:,:)
314	IF( MOD(mig(jpi-2), 3) == 0 ) p3d(jpi-1,:,:) = p3d( jpi,:,:)
315	ENDIF
316	IF( mpiRnei(nn_hls,jpso) > -1 ) THEN
317	IF( MOD(mjg( 1), 3) == 1 ) p3d(:, 1,:) = p3d(:, 2,:)
318	IF( MOD(mjg( 1), 3) == 2 ) p3d(:, 2,:) = p3d(:, 1,:)
319	ENDIF
320	IF( mpiRnei(nn_hls,jpno) > -1 ) THEN
321	IF( MOD(mjg(jpj-2), 3) == 1 ) p3d(:, jpj,:) = p3d(:,jpj-1,:)
322	IF( MOD(mjg(jpj-2), 3) == 0 ) p3d(:,jpj-1,:) = p3d(:, jpj,:)
323	ENDIF
324	CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1.0_wp )
325	ENDIF
326
327	END SUBROUTINE sum3x3_3d
328
329
330	FUNCTION glob_sum_vec_3d( cdname, ptab ) RESULT( ptmp )
331	!!----------------------------------------------------------------------
332	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
333	REAL(wp), INTENT(in) :: ptab(:,:,:) ! array on which operation is applied
334	REAL(wp), DIMENSION(SIZE(ptab,3)) :: ptmp
335	!
336	COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: ctmp
337	REAL(wp) :: ztmp
338	INTEGER :: ji , jj , jk ! dummy loop indices
339	INTEGER :: ipi, ipj, ipk ! dimensions
340	INTEGER :: iis, iie, ijs, ije ! loop start and end
341	!!-----------------------------------------------------------------------
342	!
343	ipi = SIZE(ptab,1) ! 1st dimension
344	ipj = SIZE(ptab,2) ! 2nd dimension
345	ipk = SIZE(ptab,3) ! 3rd dimension
346	!
347	IF( ipi == jpi .AND. ipj == jpj ) THEN ! do 2D loop only over the inner domain (-> avoid to use undefined values)
348	iis = Nis0 ; iie = Nie0
349	ijs = Njs0 ; ije = Nje0
350	ELSE ! I think we are never in this case...
351	iis = 1 ; iie = jpi
352	ijs = 1 ; ije = jpj
353	ENDIF
354	!
355	ALLOCATE( ctmp(ipk) )
356	!
357	DO jk = 1, ipk
358	ctmp(jk) = CMPLX( 0.e0, 0.e0, dp ) ! warning ctmp is cumulated
359	DO jj = ijs, ije
360	DO ji = iis, iie
361	ztmp = ptab(ji,jj,jk) * tmask_i(ji,jj)
362	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
363	END DO
364	END DO
365	END DO
366	CALL mpp_sum( cdname, ctmp(:) ) ! sum over the global domain
367	!
368	ptmp = REAL( ctmp(:), wp )
369	!
370	DEALLOCATE( ctmp )
371	!
372	END FUNCTION glob_sum_vec_3d
373
374	FUNCTION glob_sum_vec_4d( cdname, ptab ) RESULT( ptmp )
375	!!----------------------------------------------------------------------
376	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
377	REAL(wp), INTENT(in) :: ptab(:,:,:,:) ! array on which operation is applied
378	REAL(wp), DIMENSION(SIZE(ptab,4)) :: ptmp
379	!
380	COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: ctmp
381	REAL(wp) :: ztmp
382	INTEGER :: ji , jj , jk , jl ! dummy loop indices
383	INTEGER :: ipi, ipj, ipk, ipl ! dimensions
384	INTEGER :: iis, iie, ijs, ije ! loop start and end
385	!!-----------------------------------------------------------------------
386	!
387	ipi = SIZE(ptab,1) ! 1st dimension
388	ipj = SIZE(ptab,2) ! 2nd dimension
389	ipk = SIZE(ptab,3) ! 3rd dimension
390	ipl = SIZE(ptab,4) ! 4th dimension
391	!
392	IF( ipi == jpi .AND. ipj == jpj ) THEN ! do 2D loop only over the inner domain (-> avoid to use undefined values)
393	iis = Nis0 ; iie = Nie0
394	ijs = Njs0 ; ije = Nje0
395	ELSE ! I think we are never in this case...
396	iis = 1 ; iie = jpi
397	ijs = 1 ; ije = jpj
398	ENDIF
399	!
400	ALLOCATE( ctmp(ipl) )
401	!
402	DO jl = 1, ipl
403	ctmp(jl) = CMPLX( 0.e0, 0.e0, dp ) ! warning ctmp is cumulated
404	DO jk = 1, ipk
405	DO jj = ijs, ije
406	DO ji = iis, iie
407	ztmp = ptab(ji,jj,jk,jl) * tmask_i(ji,jj)
408	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
409	END DO
410	END DO
411	END DO
412	END DO
413	CALL mpp_sum( cdname, ctmp(:) ) ! sum over the global domain
414	!
415	ptmp = REAL( ctmp(:), wp )
416	!
417	DEALLOCATE( ctmp )
418	!
419	END FUNCTION glob_sum_vec_4d
420
421	FUNCTION glob_sum_full_vec_3d( cdname, ptab ) RESULT( ptmp )
422	!!----------------------------------------------------------------------
423	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
424	REAL(wp), INTENT(in) :: ptab(:,:,:) ! array on which operation is applied
425	REAL(wp), DIMENSION(SIZE(ptab,3)) :: ptmp
426	!
427	COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: ctmp
428	REAL(wp) :: ztmp
429	INTEGER :: ji , jj , jk ! dummy loop indices
430	INTEGER :: ipi, ipj, ipk ! dimensions
431	INTEGER :: iis, iie, ijs, ije ! loop start and end
432	!!-----------------------------------------------------------------------
433	!
434	ipi = SIZE(ptab,1) ! 1st dimension
435	ipj = SIZE(ptab,2) ! 2nd dimension
436	ipk = SIZE(ptab,3) ! 3rd dimension
437	!
438	IF( ipi == jpi .AND. ipj == jpj ) THEN ! do 2D loop only over the inner domain (-> avoid to use undefined values)
439	iis = Nis0 ; iie = Nie0
440	ijs = Njs0 ; ije = Nje0
441	ELSE ! I think we are never in this case...
442	iis = 1 ; iie = jpi
443	ijs = 1 ; ije = jpj
444	ENDIF
445	!
446	ALLOCATE( ctmp(ipk) )
447	!
448	DO jk = 1, ipk
449	ctmp(jk) = CMPLX( 0.e0, 0.e0, dp ) ! warning ctmp is cumulated
450	DO jj = ijs, ije
451	DO ji = iis, iie
452	ztmp = ptab(ji,jj,jk) * tmask_h(ji,jj)
453	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
454	END DO
455	END DO
456	END DO
457	CALL mpp_sum( cdname, ctmp(:) ) ! sum over the global domain
458	!
459	ptmp = REAL( ctmp(:), wp )
460	!
461	DEALLOCATE( ctmp )
462	!
463	END FUNCTION glob_sum_full_vec_3d
464
465	FUNCTION glob_sum_full_vec_4d( cdname, ptab ) RESULT( ptmp )
466	!!----------------------------------------------------------------------
467	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
468	REAL(wp), INTENT(in) :: ptab(:,:,:,:) ! array on which operation is applied
469	REAL(wp), DIMENSION(SIZE(ptab,4)) :: ptmp
470	!
471	COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: ctmp
472	REAL(wp) :: ztmp
473	INTEGER :: ji , jj , jk , jl ! dummy loop indices
474	INTEGER :: ipi, ipj, ipk, ipl ! dimensions
475	INTEGER :: iis, iie, ijs, ije ! loop start and end
476	!!-----------------------------------------------------------------------
477	!
478	ipi = SIZE(ptab,1) ! 1st dimension
479	ipj = SIZE(ptab,2) ! 2nd dimension
480	ipk = SIZE(ptab,3) ! 3rd dimension
481	ipl = SIZE(ptab,4) ! 4th dimension
482	!
483	IF( ipi == jpi .AND. ipj == jpj ) THEN ! do 2D loop only over the inner domain (-> avoid to use undefined values)
484	iis = Nis0 ; iie = Nie0
485	ijs = Njs0 ; ije = Nje0
486	ELSE ! I think we are never in this case...
487	iis = 1 ; iie = jpi
488	ijs = 1 ; ije = jpj
489	ENDIF
490	!
491	ALLOCATE( ctmp(ipl) )
492	!
493	DO jl = 1, ipl
494	ctmp(jl) = CMPLX( 0.e0, 0.e0, dp ) ! warning ctmp is cumulated
495	DO jk = 1, ipk
496	DO jj = ijs, ije
497	DO ji = iis, iie
498	ztmp = ptab(ji,jj,jk,jl) * tmask_h(ji,jj)
499	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
500	END DO
501	END DO
502	END DO
503	END DO
504	CALL mpp_sum( cdname, ctmp(:) ) ! sum over the global domain
505	!
506	ptmp = REAL( ctmp(:), wp )
507	!
508	DEALLOCATE( ctmp )
509	!
510	END FUNCTION glob_sum_full_vec_4d
511
512	FUNCTION glob_min_vec_3d( cdname, ptab ) RESULT( ptmp )
513	!!----------------------------------------------------------------------
514	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
515	REAL(wp), INTENT(in) :: ptab(:,:,:) ! array on which operation is applied
516	REAL(wp), DIMENSION(SIZE(ptab,3)) :: ptmp
517	!
518	INTEGER :: jk ! dummy loop indice & dimension
519	INTEGER :: ipk ! dimension
520	!!-----------------------------------------------------------------------
521	!
522	ipk = SIZE(ptab,3)
523	DO jk = 1, ipk
524	ptmp(jk) = MINVAL( ptab(:,:,jk) * tmask_i(:,:) )
525	ENDDO
526	!
527	CALL mpp_min( cdname, ptmp (:) )
528	!
529	END FUNCTION glob_min_vec_3d
530
531	FUNCTION glob_min_vec_4d( cdname, ptab ) RESULT( ptmp )
532	!!----------------------------------------------------------------------
533	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
534	REAL(wp), INTENT(in) :: ptab(:,:,:,:) ! array on which operation is applied
535	REAL(wp), DIMENSION(SIZE(ptab,4)) :: ptmp
536	!
537	INTEGER :: jk , jl ! dummy loop indice & dimension
538	INTEGER :: ipk, ipl ! dimension
539	!!-----------------------------------------------------------------------
540	!
541	ipk = SIZE(ptab,3)
542	ipl = SIZE(ptab,4)
543	DO jl = 1, ipl
544	ptmp(jl) = MINVAL( ptab(:,:,1,jl) * tmask_i(:,:) )
545	DO jk = 2, ipk
546	ptmp(jl) = MIN( ptmp(jl), MINVAL( ptab(:,:,jk,jl) * tmask_i(:,:) ) )
547	ENDDO
548	ENDDO
549	!
550	CALL mpp_min( cdname, ptmp (:) )
551	!
552	END FUNCTION glob_min_vec_4d
553
554	FUNCTION glob_max_vec_3d( cdname, ptab ) RESULT( ptmp )
555	!!----------------------------------------------------------------------
556	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
557	REAL(wp), INTENT(in) :: ptab(:,:,:) ! array on which operation is applied
558	REAL(wp), DIMENSION(SIZE(ptab,3)) :: ptmp
559	!
560	INTEGER :: jk ! dummy loop indice & dimension
561	INTEGER :: ipk ! dimension
562	!!-----------------------------------------------------------------------
563	!
564	ipk = SIZE(ptab,3)
565	DO jk = 1, ipk
566	ptmp(jk) = MAXVAL( ptab(:,:,jk) * tmask_i(:,:) )
567	ENDDO
568	!
569	CALL mpp_max( cdname, ptmp (:) )
570	!
571	END FUNCTION glob_max_vec_3d
572
573	FUNCTION glob_max_vec_4d( cdname, ptab ) RESULT( ptmp )
574	!!----------------------------------------------------------------------
575	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
576	REAL(wp), INTENT(in) :: ptab(:,:,:,:) ! array on which operation is applied
577	REAL(wp), DIMENSION(SIZE(ptab,4)) :: ptmp
578	!
579	INTEGER :: jk , jl ! dummy loop indice & dimension
580	INTEGER :: ipk, ipl ! dimension
581	!!-----------------------------------------------------------------------
582	!
583	ipk = SIZE(ptab,3)
584	ipl = SIZE(ptab,4)
585	DO jl = 1, ipl
586	ptmp(jl) = MAXVAL( ptab(:,:,1,jl) * tmask_i(:,:) )
587	DO jk = 2, ipk
588	ptmp(jl) = MAX( ptmp(jl), MAXVAL( ptab(:,:,jk,jl) * tmask_i(:,:) ) )
589	ENDDO
590	ENDDO
591	!
592	CALL mpp_max( cdname, ptmp (:) )
593	!
594	END FUNCTION glob_max_vec_4d
595
596	SUBROUTINE DDPDD( ydda, yddb )
597	!!----------------------------------------------------------------------
598	!! * ROUTINE DDPDD *
599	!!
600	!! ** Purpose : Add a scalar element to a sum
601	!!
602	!!
603	!! ** Method : The code uses the compensated summation with doublet
604	!! (sum,error) emulated useing complex numbers. ydda is the
605	!! scalar to add to the summ yddb
606	!!
607	!! ** Action : This does only work for MPI.
608	!!
609	!! References : Using Acurate Arithmetics to Improve Numerical
610	!! Reproducibility and Sability in Parallel Applications
611	!! Yun HE and Chris H. Q. DING, Journal of Supercomputing 18, 259-277, 2001
612	!!----------------------------------------------------------------------
613	COMPLEX(dp), INTENT(in ) :: ydda
614	COMPLEX(dp), INTENT(inout) :: yddb
615	!
616	REAL(dp) :: zerr, zt1, zt2 ! local work variables
617	!!-----------------------------------------------------------------------
618	!
619	! Compute ydda + yddb using Knuth's trick.
620	zt1 = REAL(ydda) + REAL(yddb)
621	zerr = zt1 - REAL(ydda)
622	zt2 = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) ) &
623	& + AIMAG(ydda) + AIMAG(yddb)
624	!
625	! The result is t1 + t2, after normalization.
626	yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
627	!
628	END SUBROUTINE DDPDD
629
630	#if defined key_nosignedzero
631	!!----------------------------------------------------------------------
632	!! 'key_nosignedzero' F90 SIGN
633	!!----------------------------------------------------------------------
634
635	FUNCTION SIGN_SCALAR( pa, pb )
636	!!-----------------------------------------------------------------------
637	!! * FUNCTION SIGN_SCALAR *
638	!!
639	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
640	!!-----------------------------------------------------------------------
641	REAL(wp) :: pa,pb ! input
642	REAL(wp) :: SIGN_SCALAR ! result
643	!!-----------------------------------------------------------------------
644	IF ( pb >= 0.e0) THEN ; SIGN_SCALAR = ABS(pa)
645	ELSE ; SIGN_SCALAR =-ABS(pa)
646	ENDIF
647	END FUNCTION SIGN_SCALAR
648
649
650	FUNCTION SIGN_ARRAY_1D( pa, pb )
651	!!-----------------------------------------------------------------------
652	!! * FUNCTION SIGN_ARRAY_1D *
653	!!
654	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
655	!!-----------------------------------------------------------------------
656	REAL(wp) :: pa,pb(:) ! input
657	REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1)) ! result
658	!!-----------------------------------------------------------------------
659	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_1D = ABS(pa)
660	ELSEWHERE ; SIGN_ARRAY_1D =-ABS(pa)
661	END WHERE
662	END FUNCTION SIGN_ARRAY_1D
663
664
665	FUNCTION SIGN_ARRAY_2D(pa,pb)
666	!!-----------------------------------------------------------------------
667	!! * FUNCTION SIGN_ARRAY_2D *
668	!!
669	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
670	!!-----------------------------------------------------------------------
671	REAL(wp) :: pa,pb(:,:) ! input
672	REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2)) ! result
673	!!-----------------------------------------------------------------------
674	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_2D = ABS(pa)
675	ELSEWHERE ; SIGN_ARRAY_2D =-ABS(pa)
676	END WHERE
677	END FUNCTION SIGN_ARRAY_2D
678
679	FUNCTION SIGN_ARRAY_3D(pa,pb)
680	!!-----------------------------------------------------------------------
681	!! * FUNCTION SIGN_ARRAY_3D *
682	!!
683	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
684	!!-----------------------------------------------------------------------
685	REAL(wp) :: pa,pb(:,:,:) ! input
686	REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
687	!!-----------------------------------------------------------------------
688	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_3D = ABS(pa)
689	ELSEWHERE ; SIGN_ARRAY_3D =-ABS(pa)
690	END WHERE
691	END FUNCTION SIGN_ARRAY_3D
692
693
694	FUNCTION SIGN_ARRAY_1D_A(pa,pb)
695	!!-----------------------------------------------------------------------
696	!! * FUNCTION SIGN_ARRAY_1D_A *
697	!!
698	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
699	!!-----------------------------------------------------------------------
700	REAL(wp) :: pa(:),pb(:) ! input
701	REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(pb,1)) ! result
702	!!-----------------------------------------------------------------------
703	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_1D_A = ABS(pa)
704	ELSEWHERE ; SIGN_ARRAY_1D_A =-ABS(pa)
705	END WHERE
706	END FUNCTION SIGN_ARRAY_1D_A
707
708
709	FUNCTION SIGN_ARRAY_2D_A(pa,pb)
710	!!-----------------------------------------------------------------------
711	!! * FUNCTION SIGN_ARRAY_2D_A *
712	!!
713	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
714	!!-----------------------------------------------------------------------
715	REAL(wp) :: pa(:,:),pb(:,:) ! input
716	REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2)) ! result
717	!!-----------------------------------------------------------------------
718	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_2D_A = ABS(pa)
719	ELSEWHERE ; SIGN_ARRAY_2D_A =-ABS(pa)
720	END WHERE
721	END FUNCTION SIGN_ARRAY_2D_A
722
723
724	FUNCTION SIGN_ARRAY_3D_A(pa,pb)
725	!!-----------------------------------------------------------------------
726	!! * FUNCTION SIGN_ARRAY_3D_A *
727	!!
728	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
729	!!-----------------------------------------------------------------------
730	REAL(wp) :: pa(:,:,:),pb(:,:,:) ! input
731	REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
732	!!-----------------------------------------------------------------------
733	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_3D_A = ABS(pa)
734	ELSEWHERE ; SIGN_ARRAY_3D_A =-ABS(pa)
735	END WHERE
736	END FUNCTION SIGN_ARRAY_3D_A
737
738
739	FUNCTION SIGN_ARRAY_1D_B(pa,pb)
740	!!-----------------------------------------------------------------------
741	!! * FUNCTION SIGN_ARRAY_1D_B *
742	!!
743	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
744	!!-----------------------------------------------------------------------
745	REAL(wp) :: pa(:),pb ! input
746	REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1)) ! result
747	!!-----------------------------------------------------------------------
748	IF( pb >= 0.e0 ) THEN ; SIGN_ARRAY_1D_B = ABS(pa)
749	ELSE ; SIGN_ARRAY_1D_B =-ABS(pa)
750	ENDIF
751	END FUNCTION SIGN_ARRAY_1D_B
752
753
754	FUNCTION SIGN_ARRAY_2D_B(pa,pb)
755	!!-----------------------------------------------------------------------
756	!! * FUNCTION SIGN_ARRAY_2D_B *
757	!!
758	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
759	!!-----------------------------------------------------------------------
760	REAL(wp) :: pa(:,:),pb ! input
761	REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2)) ! result
762	!!-----------------------------------------------------------------------
763	IF( pb >= 0.e0 ) THEN ; SIGN_ARRAY_2D_B = ABS(pa)
764	ELSE ; SIGN_ARRAY_2D_B =-ABS(pa)
765	ENDIF
766	END FUNCTION SIGN_ARRAY_2D_B
767
768
769	FUNCTION SIGN_ARRAY_3D_B(pa,pb)
770	!!-----------------------------------------------------------------------
771	!! * FUNCTION SIGN_ARRAY_3D_B *
772	!!
773	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
774	!!-----------------------------------------------------------------------
775	REAL(wp) :: pa(:,:,:),pb ! input
776	REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3)) ! result
777	!!-----------------------------------------------------------------------
778	IF( pb >= 0.e0 ) THEN ; SIGN_ARRAY_3D_B = ABS(pa)
779	ELSE ; SIGN_ARRAY_3D_B =-ABS(pa)
780	ENDIF
781	END FUNCTION SIGN_ARRAY_3D_B
782	#endif
783
784	!!======================================================================
785	END MODULE lib_fortran

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

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

Context Navigation

source: NEMO/trunk/src/OCE/lib_fortran.F90

Download in other formats: