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.

lib_fortran.F90 in NEMO/trunk/src/OCE – NEMO

Context Navigation

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

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

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

Download in other formats: