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/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/OCE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/OCE/lib_fortran.F90 @ 13307

Last change on this file since 13307 was 13090, checked in by andmirek, 4 years ago

Ticket #2482: calculate sums in the same way as in #2197, commit 11069

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