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

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

Last change on this file since 15376 was 15376, checked in by clem, 12 months 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 
1MODULE 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   !!----------------------------------------------------------------------
87CONTAINS
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   !!======================================================================
785END MODULE lib_fortran
Note: See TracBrowser for help on using the repository browser.