source: branches/DEV_1879_mpp_rep/NEMO/OPA_SRC/lib_fortran.F90 @ 2003

Last change on this file since 2003 was 2003, checked in by rblod, 10 years ago

Add missing module for mpp reproducibility, see ticket #677

File size: 22.9 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  Michael Dunphy, Rachid BENSHILA Original code
7   !!----------------------------------------------------------------------
8   !!   OPA 9.0 , LOCEAN-IPSL (2005)
9   !! $Id: $
10   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
11   !!----------------------------------------------------------------------
12   USE par_oce 
13   USE par_kind
14   USE lib_mpp         ! distributed memory computing
15   USE dom_oce 
16   USE in_out_manager
17
18   IMPLICIT NONE
19   PRIVATE
20
21   PUBLIC glob_sum
22#if defined key_nosignedzeo
23   PUBLIC SIGN
24#endif
25
26   INTERFACE glob_sum
27#if defined key_mpp_rep1
28      MODULE PROCEDURE mpp_sum_indep
29#elif defined key_mpp_rep2     
30      MODULE PROCEDURE mpp_sum_cmpx
31#else
32      MODULE PROCEDURE glob_sum_2d, glob_sum_3d,glob_sum_2d_a, glob_sum_3d_a 
33#endif
34   END INTERFACE
35
36#if defined key_nosignedzeo   
37   INTERFACE SIGN
38      MODULE PROCEDURE SIGN_SCALAR, SIGN_ARRAY_1D, SIGN_ARRAY_2D, SIGN_ARRAY_3D, &
39                       SIGN_ARRAY_1D_A, SIGN_ARRAY_2D_A, SIGN_ARRAY_3D_A,  & 
40                       SIGN_ARRAY_1D_B, SIGN_ARRAY_2D_B, SIGN_ARRAY_3D_B 
41   END INTERFACE
42#endif
43
44CONTAINS
45
46   FUNCTION glob_sum_2d( ptab )
47      !!-----------------------------------------------------------------------
48      !!                  ***  FUNCTION  glob_sum_2D  ***
49      !!
50      !! ** Purpose : perform a sum on the global domain of a 2D array
51      !!-----------------------------------------------------------------------
52      REAL(wp), DIMENSION(:,:),INTENT(in) :: ptab
53      REAL(wp) :: glob_sum_2d
54      !!-----------------------------------------------------------------------
55
56      glob_sum_2d = SUM( ptab(:,:)*tmask_i(:,:) )
57      IF( lk_mpp )   CALL mpp_sum( glob_sum_2d )
58       
59   END FUNCTION glob_sum_2d
60 
61   FUNCTION glob_sum_3d( ptab )
62      !!-----------------------------------------------------------------------
63      !!                  ***  FUNCTION  glob_sum_3D  ***
64      !!
65      !! ** Purpose : perform a sum on the global domain of a 3D array
66      !!-----------------------------------------------------------------------
67      REAL(wp), DIMENSION(:,:,:) :: ptab
68      REAL(wp) :: glob_sum_3d
69      !
70      INTEGER :: jk
71      !!-----------------------------------------------------------------------
72       
73      GLOB_SUM_3D = 0.e0
74      DO jk = 1, jpk
75         glob_sum_3d = glob_sum_3d + SUM( ptab(:,:,jk)*tmask_i(:,:) )
76      END DO
77      IF( lk_mpp )   CALL mpp_sum( glob_sum_3d )
78       
79   END FUNCTION glob_sum_3d
80
81   FUNCTION glob_sum_2d_a( ptab1, ptab2 )
82      !!-----------------------------------------------------------------------
83      !!                  ***  FUNCTION  glob_sum_2D _a ***
84      !!
85      !! ** Purpose : perform a sum on the global domain of two 2D array
86      !!-----------------------------------------------------------------------
87      REAL(wp), DIMENSION(:,:) :: ptab1, ptab2
88      REAL(wp), DIMENSION(2)   :: glob_sum_2d_a
89      !!-----------------------------------------------------------------------
90                   
91      glob_sum_2d_a(1) = SUM( ptab1(:,:)*tmask_i(:,:) )
92      glob_sum_2d_a(2) = SUM( ptab2(:,:)*tmask_i(:,:) )
93      IF( lk_mpp )   CALL mpp_sum( glob_sum_2d_a,2 )
94       
95   END FUNCTION glob_sum_2d_a
96 
97   FUNCTION glob_sum_3d_a( ptab1, ptab2 )
98      !!-----------------------------------------------------------------------
99      !!                  ***  FUNCTION  glob_sum_3D_a ***
100      !!
101      !! ** Purpose : perform a sum on the global domain of two 3D array
102      !!-----------------------------------------------------------------------
103      REAL(wp), DIMENSION(:,:,:) :: ptab1, ptab2
104      REAL(wp), DIMENSION(2)     :: glob_sum_3d_a
105      !
106      INTEGER :: jk
107      !!-----------------------------------------------------------------------
108       
109      glob_sum_3d_a(:) = 0.e0
110      DO jk = 1, jpk
111         glob_sum_3d_a(1) = glob_sum_3d_a(1) + SUM( ptab1(:,:,jk)*tmask_i(:,:) )
112         glob_sum_3d_a(2) = glob_sum_3d_a(2) + SUM( ptab2(:,:,jk)*tmask_i(:,:) )
113      END DO
114      IF( lk_mpp )   CALL mpp_sum( glob_sum_3d_a,2 )
115       
116   END FUNCTION glob_sum_3d_a
117
118#if defined key_mpp_rep2 
119   FUNCTION mpp_sum_cmpx( pval )
120      !!----------------------------------------------------------------------
121      !!                  ***  FUNCTION  mpp_sum_cmpx ***
122      !!
123      !! ** Purpose : perform a sum in calling DDPDD routine
124      !!
125      !!----------------------------------------------------------------------
126      REAL(wp) :: mpp_sum_cmpx
127      !
128      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: &
129         & pval
130      COMPLEX(wp):: ctmp
131      REAL(wp) ::ztmp
132      INTEGER :: ji,jj
133      !!-----------------------------------------------------------------------
134     
135      ztmp = 0.e0
136      ctmp = CMPLX(0.e0,0.e0,wp)
137      DO jj = 1,jpj
138         DO ji =1, jpi
139         ztmp =  pval(ji,jj) * tmask_i(ji,jj)
140         CALL DDPDD(CMPLX(ztmp,0.e0,wp),ctmp)
141         END DO
142      END DO
143      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
144      mpp_sum_cmpx= REAL(ctmp,wp)
145       
146   END FUNCTION mpp_sum_cmpx
147
148   SUBROUTINE DDPDD( ydda, yddb )
149      !!----------------------------------------------------------------------
150      !!               ***  ROUTINE DDPDD ***
151      !!         
152      !! ** Purpose : Add a scalar element to a sum
153      !!             
154      !!
155      !! ** Method  : The code uses the compensated summation with doublet
156      !!              (sum,error) emulated useing complex numbers. ydda is the
157      !!               scalar to add to the summ yddb
158      !!
159      !! ** Action  : This does only work for MPI.
160      !!
161      !! References : Using Acurate Arithmetics to Improve Numerical
162      !!              Reproducibility and Sability in Parallel Applications
163      !!              Yun HE and Chris H. Q. DING, Journal of Supercomputing
164      !!                                            18, 259-277, 2001
165      !!----------------------------------------------------------------------
166
167      COMPLEX(wp), INTENT(in)     :: ydda
168      COMPLEX(wp), INTENT(inout)  :: yddb
169     
170      REAL(wp) :: zerr, zt1, zt2  ! local work variables
171
172      ! Compute ydda + yddb using Knuth's trick.
173      zt1  = real(ydda) + real(yddb)
174      zerr = zt1 - real(ydda)
175      zt2  = ((real(yddb) - zerr) + (real(ydda) - (zt1 - zerr))) &
176            + aimag(ydda) + aimag(yddb)
177
178      ! The result is t1 + t2, after normalization.
179      yddb = cmplx ( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1),wp )
180     
181   END SUBROUTINE DDPDD
182#endif
183
184#if defined key_mpp_rep1
185   FUNCTION mpp_sum_indep( pval )
186      !!----------------------------------------------------------------------
187      !!               ***  ROUTINE mpp_sum_indep ***
188      !!         
189      !! ** Purpose : Sum all elements in the pval array in
190      !!              an accurate order-independent way.
191      !!
192      !! ** Method  : The code iterates the compensated summation until the
193      !!              result is guaranteed to be within 4*eps of the true sum.
194      !!              It then rounds the result to the nearest floating-point
195      !!              number whose last three bits are zero, thereby
196      !!              guaranteeing an order-independent result.
197      !!
198      !! ** Action  : This does only work for MPI.
199      !!              It does not work for SHMEM.      !!
200      !! References : M. Fisher (ECMWF): IFS code + personal communication
201      !!              The algorithm is based on Ogita et al. (2005)
202      !!              SIAM J. Sci. Computing, Vol.26, No.6, pp1955-1988.
203      !!              This is based in turn on an algorithm
204      !!              by Knuth (1969, seminumerical algorithms).
205      !!
206      !! History :
207      !!        !  07-07  (K. Mogensen)  Original code heavily based on IFS.
208      !!----------------------------------------------------------------------
209      REAL(wp) mpp_sum_indep
210      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pval
211     
212      REAL(wp), DIMENSION(3) :: zbuffl
213      REAL(wp), DIMENSION(:), ALLOCATABLE :: zpsums, zperrs, zpcors, zbuffg, zp
214      REAL(wp) :: zcorr, zerr,  zolderr, zbeta, zres
215      INTEGER, DIMENSION(:), allocatable :: irecv, istart
216      INTEGER :: ikn, jj
217
218      ! initialise to avoid uninitialised variables trapping of some compilers to complain.
219      zres = 0.0_wp ; zerr = 0.0_wp ; zbuffl(:) = 0.0_wp
220      ! Get global number of elements
221      ikn = SIZE(pval)
222# ifdef key_mpp
223      CALL mpp_sum( ikn )
224# endif
225      ! Check that the the algorithm can work
226     
227      IF ( ( REAL( 2 * ikn ) * EPSILON( zres ) ) >= 1.0 ) THEN
228         CALL ctl_stop('mpp_sum_indep:', &
229            &          'size of array is too large to guarantee error bounds')
230      ENDIF
231     
232      ALLOCATE( &
233         & zp(MAX(ikn,1)),                & 
234         & zbuffg(jpnij*SIZE(zbuffl)), &
235         & zpsums(jpnij),              &
236         & zperrs(jpnij),              &
237         & zpcors(jpnij)               &
238         & )
239
240      zolderr = HUGE(zerr)
241
242      ! Copy the input array. This avoids some tricky indexing, at the
243      ! expense of some inefficency.
244
245      IF ( ikn > 0 ) THEN
246         zp(:) = RESHAPE(pval, (/ jpi * jpj /) )
247      ELSE
248         zp(1) = 0.0_wp
249      ENDIF
250     
251      k_loop: DO
252
253         ! Transform local arrays
254         
255         IF ( ikn > 0 ) THEN
256            CALL comp_sum ( zp, ikn, zcorr, zerr )
257         ENDIF
258
259         ! Gather partial sums and error bounds to all processors
260
261         zbuffl(1) = zp(MAX(ikn,1))
262         
263         IF ( ikn > 0 ) THEN
264            zbuffl(2) = zerr
265            zbuffl(3) = zcorr
266         ELSE           
267            zbuffl(2) = 0.0_wp
268            zbuffl(3) = 0.0_wp
269         ENDIF
270         
271         IF ( jpnij > 1 ) THEN
272            ALLOCATE( &
273               & irecv(jpnij), &
274               & istart(jpnij) &
275               & )
276            CALL mpp_allgatherv( zbuffl, SIZE(zbuffl), &
277               & zbuffg, jpnij * SIZE(zbuffl), irecv, istart )
278            DEALLOCATE( &
279               & irecv, &
280               & istart &
281               & )
282           
283            DO jj = 1, jpnij
284               zpsums(jj) = zbuffg(1+(jj-1)*SIZE(zbuffl))
285               zperrs(jj) = zbuffg(2+(jj-1)*SIZE(zbuffl))
286               zpcors(jj) = zbuffg(3+(jj-1)*SIZE(zbuffl))
287            END DO
288
289         ELSE
290            zpsums(1) = zbuffl(1)
291            zperrs(1) = zbuffl(2)
292            zpcors(1) = zbuffl(3)
293         ENDIF
294
295         ! Transform partial sums
296         CALL comp_sum( zpsums, jpnij, zcorr, zerr )
297         zerr  = zerr  + SUM(zperrs)
298         zcorr = zcorr + SUM(zpcors)
299         
300         ! Calculate final result         
301         zres = zpsums(jpnij) + zcorr
302
303         ! Calculate error bound. This is corollary 4.7 from Ogita et al.
304         ! (2005)         
305         zbeta = zerr *(  REAL( 2*ikn, wp ) * EPSILON(zres) ) &
306            &  /(1.0_wp - REAL( 2*ikn, wp ) * EPSILON(zres) )
307
308         zerr = EPSILON(zres) * ABS(zres) &
309            & +(zbeta + ( 2.0_wp * EPSILON(zres) * EPSILON(zres) * ABS(zres) &
310            &            +3.0_wp * TINY(zres) ) )
311
312         ! Update the last element of the local array
313         zp(MAX(ikn,1)) = zpsums(nproc+1)
314
315         ! Exit if the global error is small enough
316         IF ( zerr < 4.0_wp * SPACING(zres) ) EXIT k_loop
317         
318         ! Take appropriate action if ZRES cannot be sufficiently refined.         
319         IF (zerr >= zolderr) THEN
320            CALL ctl_stop('Failed to refine sum', &
321               &          'Warning: Possiblity of non-reproducible results')
322         ENDIF
323
324         zolderr = zerr
325         
326      ENDDO k_loop
327
328      ! At this stage, we have guaranteed that ZRES less than 4*EPS
329      ! away from the exact sum. There are only four floating point
330      ! numbers in this range. So, if we find the nearest number that
331      ! has its last three bits zero, then we have a reproducible result.
332
333      mpp_sum_indep = fround(zres)
334
335      DEALLOCATE( &
336         & zpcors, &
337         & zperrs, &
338         & zpsums, &
339         & zbuffg, &
340         & zp      &
341         & )
342
343   END FUNCTION mpp_sum_indep
344
345   SUBROUTINE comp_sum( pval, kn, pcorr, perr )
346      !!----------------------------------------------------------------------
347      !!               ***  ROUTINE comp_sum ***
348      !!         
349      !! ** Purpose : To perform compensated (i.e. accurate) summation.
350      !!
351      !! ** Method  : These routines transform the elements of the array P,
352      !!              such that:
353      !!              1)  pval(kn)         contains sum(pval)
354      !!              2)  pval(1)...pval(kn-1) contain the rounding errors
355      !!                  that were made in calculating sum(pval).
356      !!              3)  The exact sum of the elements of pval is unmodified.
357      !!              On return, pcorr contains the sum of the rounding errors,
358      !!              perr contains the sum of their absolute values.
359      !!              After calling this routine, an accurate sum of the
360      !!              elements of pval can be calculated as res=pval(n)+pcorr.
361      !!
362      !! ** Action  :
363      !!
364      !! References : M. Fisher (ECMWF) IFS code + personal communications
365      !!
366      !! History :
367      !!        !  07-07  (K. Mogensen)  Original code heavily based on IFS
368      !!--------------------------------------------------------------------
369      INTEGER, INTENT(IN) :: kn         ! Number of elements in input array
370      REAL(wp), DIMENSION(kn), INTENT(INOUT) :: pval    ! Input array to be sum on input
371                                                        ! pval(kn) = sum (pval) on output
372                                                        ! pval(1)...pval(kn-1) = rounding errors on output
373      REAL(wp) :: pcorr   ! Sum of rounding errors
374      REAL(wp) :: perr       ! Sum of absolute rounding errors
375      !! * Local declarations
376      REAL(wp) :: zx, zz, zpsum
377      INTEGER :: jj
378
379      pcorr = 0.0_wp
380      perr  = 0.0_wp
381     
382      zpsum = pval(1)
383     
384      DO jj = 2, kn
385
386         ! It is vital that these 4 lines are not optimized in any way that
387         ! changes the results.
388         zx         = pval(jj) + zpsum
389         zz         = zx - pval(jj)
390         pval(jj-1) = ( pval(jj) - ( zx - zz ) ) + ( zpsum - zz ) 
391         zpsum      = zx
392
393         ! Accumulate the correction and the error         
394         pcorr      = pcorr + pval(jj-1)
395         perr       = perr + ABS( pval(jj-1) )
396
397      END DO
398
399      pval(kn) = zpsum
400     
401   END SUBROUTINE comp_sum
402
403   FUNCTION fround(pres)
404      !!----------------------------------------------------------------------
405      !!               ***  ROUTINE fround ***
406      !!         
407      !! ** Purpose : Rounding of floating-point number
408      !!
409      !! ** Method  : Returns the value of PRES rounded to the nearest
410      !!              floating-point number that has its last three bits zero
411      !!              This works on big-endian and little-endian machines.
412      !!
413      !! ** Action  :
414      !!
415      !! References : M. Fisher (ECMWF) IFS code + personal communication
416      !!
417      !! History :
418      !!        !  07-07  (K. Mogensen)  Original code heavily based on IFS.
419      !!----------------------------------------------------------------------
420      REAL(wp) fround
421      REAL(wp), INTENT(IN) :: pres      ! Value to be rounded
422      !
423      REAL(wp) :: zz(2), zup, zdown
424      INTEGER  :: ii(2), iequiv(8), ints_per_real, i_low_word
425      INTEGER  :: jj
426
427      ii(:) = 1
428      zz(:) = 1.0_wp
429
430      ! Warning: If wp = 64 bits (or 32 bits for key_sp) this will not work.
431
432#if defined key_sp
433      ints_per_real = 32 / BIT_SIZE(ii)
434#else
435      ints_per_real = 64 / BIT_SIZE(ii)
436#endif
437
438      ! Test whether big-endian or little-endian
439     
440      zup = -1.0_wp
441      iequiv(1:ints_per_real) = TRANSFER(zup,iequiv(1:ints_per_real))
442     
443      IF ( iequiv(1) == 0 ) THEN
444         i_low_word = 1                ! Little-endian
445      ELSE
446         i_low_word = ints_per_real    ! Big-endian
447      ENDIF
448
449      ! Find the nearest number with all 3 lowest-order bits zeroed
450     
451      iequiv(1:ints_per_real) = transfer(pres,iequiv(1:ints_per_real))
452      zup    = pres
453      zdown  = pres
454
455      IF (IBITS(iequiv(i_low_word),0,3)/=0) THEN
456
457         DO jj = 1, 4
458
459            zup = NEAREST( zup, 1.0_wp )
460            iequiv(1:ints_per_real) = TRANSFER( zup, iequiv(1:ints_per_real) )
461            IF ( IBITS( iequiv(i_low_word), 0, 3 ) == 0 ) EXIT
462
463            zdown = NEAREST( zdown, -1.0 )
464            iequiv(1:ints_per_real) = TRANSFER( zdown, iequiv(1:ints_per_real))
465            IF ( IBITS( iequiv(i_low_word),0,3) == 0 ) EXIT
466
467         END DO
468
469         IF ( IBITS( iequiv( i_low_word ), 0, 3) /= 0 ) THEN
470            CALL ctl_stop('Fround:','This is not possible')
471         ENDIF
472         
473      ENDIF
474     
475      fround = TRANSFER( iequiv(1:ints_per_real), pres )
476     
477   END FUNCTION fround
478#endif
479
480
481#if defined key_nosignedzero
482   FUNCTION SIGN_SCALAR(pa,pb)
483      !!-----------------------------------------------------------------------
484      !!                  ***  FUNCTION SIGN_SCALAR  ***
485      !!
486      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
487      !!-----------------------------------------------------------------------
488      REAL(wp) :: pa,pb          ! input
489      REAL(wp) :: SIGN_SCALAR  ! result
490      IF ( pb >= 0.e0) THEN
491         SIGN_SCALAR = ABS(pa)
492      ELSE
493         SIGN_SCALAR =-ABS(pa)
494      ENDIF
495
496   END FUNCTION SIGN_SCALAR
497
498   FUNCTION SIGN_ARRAY_1D(pa,pb) 
499      !!-----------------------------------------------------------------------
500      !!                  ***  FUNCTION SIGN_ARRAY_1D  ***
501      !!
502      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
503      !!-----------------------------------------------------------------------
504      REAL(wp) :: pa,pb(:)      ! input
505      REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1))  ! result
506      WHERE ( pb >= 0.e0 )
507         SIGN_ARRAY_1D = ABS(pa)
508      ELSEWHERE
509         SIGN_ARRAY_1D =-ABS(pa)
510      END WHERE
511
512   END FUNCTION SIGN_ARRAY_1D
513
514   FUNCTION SIGN_ARRAY_2D(pa,pb) 
515      !!-----------------------------------------------------------------------
516      !!                  ***  FUNCTION SIGN_ARRAY_2D  ***
517      !!
518      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
519      !!-----------------------------------------------------------------------
520      REAL(wp) :: pa,pb(:,:)      ! input
521      REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2))  ! result
522
523      WHERE ( pb >= 0.e0 )
524         SIGN_ARRAY_2D = ABS(pa)
525      ELSEWHERE
526         SIGN_ARRAY_2D =-ABS(pa)
527      END WHERE
528
529   END FUNCTION SIGN_ARRAY_2D
530
531   FUNCTION SIGN_ARRAY_3D(pa,pb) 
532      !!-----------------------------------------------------------------------
533      !!                  ***  FUNCTION SIGN_ARRAY_3D  ***
534      !!
535      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
536      !!-----------------------------------------------------------------------
537      REAL(wp) :: pa,pb(:,:,:)      ! input
538      REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3))  ! result
539      WHERE ( pb >= 0.e0 )
540         SIGN_ARRAY_3D = ABS(pa)
541      ELSEWHERE
542         SIGN_ARRAY_3D =-ABS(pa)
543      END WHERE
544
545   END FUNCTION SIGN_ARRAY_3D
546
547   FUNCTION SIGN_ARRAY_1D_A(pa,pb) 
548      !!-----------------------------------------------------------------------
549      !!                  ***  FUNCTION SIGN_ARRAY_1D_A  ***
550      !!
551      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
552      !!-----------------------------------------------------------------------
553      REAL(wp) :: pa(:),pb(:)      ! input
554      REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(b,1))  ! result
555
556      WHERE ( pb >= 0.e0 )
557         SIGN_ARRAY_1D_A = ABS(pa)
558      ELSEWHERE
559         SIGN_ARRAY_1D_A =-ABS(pa)
560      END WHERE
561
562   END FUNCTION SIGN_ARRAY_1D_A
563
564   FUNCTION SIGN_ARRAY_2D_A(pa,pb) 
565      !!-----------------------------------------------------------------------
566      !!                  ***  FUNCTION SIGN_ARRAY_2D_A  ***
567      !!
568      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
569      !!-----------------------------------------------------------------------
570      REAL(wp) :: pa(:,:),pb(:,:)      ! input
571      REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2))  ! result
572
573      WHERE ( pb >= 0.e0 )
574         SIGN_ARRAY_2D_A = ABS(pa)
575      ELSEWHERE
576         SIGN_ARRAY_2D_A =-ABS(pa)
577      END WHERE
578
579   END FUNCTION SIGN_ARRAY_2D_A
580
581   FUNCTION SIGN_ARRAY_3D_A(pa,pb) 
582      !!-----------------------------------------------------------------------
583      !!                  ***  FUNCTION SIGN_ARRAY_3D_A  ***
584      !!
585      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
586      !!-----------------------------------------------------------------------
587      REAL(wp) :: pa(:,:,:),pb(:,:,:)  ! input
588      REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
589
590      WHERE ( pb >= 0.e0 )
591         SIGN_ARRAY_3D_A = ABS(pa)
592      ELSEWHERE
593         SIGN_ARRAY_3D_A =-ABS(pa)
594      END WHERE
595
596   END FUNCTION SIGN_ARRAY_3D_A
597
598   FUNCTION SIGN_ARRAY_1D_B(pa,pb) 
599      !!-----------------------------------------------------------------------
600      !!                  ***  FUNCTION SIGN_ARRAY_1D_B  ***
601      !!
602      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
603      !!-----------------------------------------------------------------------
604      REAL(wp) :: pa(:),pb      ! input
605      REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1))  ! result
606
607      IF ( pb >= 0.e0 ) THEN
608         SIGN_ARRAY_1D_B = ABS(pa)
609      ELSE
610         SIGN_ARRAY_1D_B =-ABS(pa)
611      ENDIF
612
613   END FUNCTION SIGN_ARRAY_1D_B
614
615   FUNCTION SIGN_ARRAY_2D_B(pa,pb) 
616      !!-----------------------------------------------------------------------
617      !!                  ***  FUNCTION SIGN_ARRAY_2D_B  ***
618      !!
619      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
620      !!-----------------------------------------------------------------------
621      REAL(wp) :: pa(:,:),pb      ! input
622      REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2))  ! result
623
624      IF ( pb >= 0.e0 ) THEN
625         SIGN_ARRAY_2D_B = ABS(pa)
626      ELSE
627         SIGN_ARRAY_2D_B =-ABS(pa)
628      ENDIF
629
630   END FUNCTION SIGN_ARRAY_2D_B
631
632   FUNCTION SIGN_ARRAY_3D_B(pa,pb) 
633      !!-----------------------------------------------------------------------
634      !!                  ***  FUNCTION SIGN_ARRAY_3D_B  ***
635      !!
636      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
637      !!-----------------------------------------------------------------------
638      REAL(wp) :: pa(:,:,:),pb      ! input
639      REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3))  ! result
640
641      IF (pb >= 0.e0 ) THEN
642         SIGN_ARRAY_3D_B = ABS(pa)
643      ELSE
644         SIGN_ARRAY_3D_B =-ABS(pa)
645      ENDIF
646 
647   END FUNCTION SIGN_ARRAY_3D_B
648#endif
649
650END MODULE lib_fortran
Note: See TracBrowser for help on using the repository browser.