ungrib/src/ngl/g2/intmath.f

   1       module intmath
   2       implicit none
   3
   4       interface ilog2
   5       ! log(x)/log(2)
   6       module procedure ilog2_8
   7       module procedure ilog2_4
   8       module procedure ilog2_2
   9       module procedure ilog2_1
  10       end interface ilog2
  11
  12       interface i1log2
  13       ! log(x+1)/log(2) unless x=maxint, in which case log(x)/log(2)
  14       module procedure i1log2_8
  15       module procedure i1log2_4
  16       module procedure i1log2_2
  17       module procedure i1log2_1
  18       end interface i1log2
  19
  20       contains
  21
  22       ! ----------------------------------------------------------------
  23
  24       function i1log2_8(ival)
  25       implicit none
  26       integer(kind=8), value :: ival
  27       integer(kind=8)::i1log2_8
  28       integer(kind=8), parameter :: one=1
  29       if(ival+one<ival) then
  30          i1log2_8=ilog2_8(ival)
  31       else
  32          i1log2_8=ilog2_8(ival+one)
  33       endif
  34       end function i1log2_8
  35
  36       ! ----------------------------------------------------------------
  37
  38       function i1log2_4(ival)
  39       implicit none
  40       integer(kind=4), value :: ival
  41       integer(kind=4)::i1log2_4
  42       integer(kind=4), parameter :: one=1
  43       if(ival+one<ival) then
  44          i1log2_4=ilog2_4(ival)
  45       else
  46          i1log2_4=ilog2_4(ival+one)
  47       endif
  48       end function i1log2_4
  49
  50       ! ----------------------------------------------------------------
  51
  52       function i1log2_2(ival)
  53       implicit none
  54       integer(kind=2), value :: ival
  55       integer(kind=2)::i1log2_2
  56       integer(kind=2), parameter :: one=1
  57       if(ival+one<ival) then
  58          i1log2_2=ilog2_2(ival)
  59       else
  60          i1log2_2=ilog2_2(ival+one)
  61       endif
  62       end function i1log2_2
  63
  64       ! ----------------------------------------------------------------
  65
  66       function i1log2_1(ival)
  67       implicit none
  68       integer(kind=1), value :: ival
  69       integer(kind=1)::i1log2_1
  70       integer(kind=1), parameter :: one=1
  71       if(ival+one<ival) then
  72          i1log2_1=ilog2_1(ival)
  73       else
  74          i1log2_1=ilog2_1(ival+one)
  75       endif
  76       end function i1log2_1
  77
  78       ! ----------------------------------------------------------------
  79
  80       function ilog2_8(i_in)
  81       implicit none
  82       integer(kind=8), value :: i_in
  83       integer(kind=8)::ilog2_8,i
  84       ilog2_8=0
  85       i=i_in
  86       if(i<=0) return
  87       if(iand(i,i-1)/=0) then
  88          !write(0,*) 'iand i-1'
  89          ilog2_8=1
  90       endif
  91       if(iand(i,Z'FFFFFFFF00000000')/=0) then
  92          ilog2_8=ilog2_8+32
  93          i=ishft(i,-32)
  94          !write(0,*) 'iand ffffffff',i,ilog2_8
  95       endif
  96       if(iand(i,Z'00000000FFFF0000')/=0) then
  97          ilog2_8=ilog2_8+16
  98          i=ishft(i,-16)
  99          !write(0,*) 'iand ffff' ,i,ilog2_8
 100       endif
 101       if(iand(i,Z'000000000000FF00')/=0) then
 102          ilog2_8=ilog2_8+8
 103          i=ishft(i,-8)
 104          !write(0,*) 'iand ff',i,ilog2_8
 105       endif
 106       if(iand(i,Z'00000000000000F0')/=0) then
 107          ilog2_8=ilog2_8+4
 108          i=ishft(i,-4)
 109          !write(0,*) 'iand f',i,ilog2_8
 110       endif
 111       if(iand(i,Z'000000000000000C')/=0) then
 112          ilog2_8=ilog2_8+2
 113          i=ishft(i,-2)
 114          !write(0,*) 'iand c',i,ilog2_8
 115       endif
 116       if(iand(i,Z'0000000000000002')/=0) then
 117          ilog2_8=ilog2_8+1
 118          i=ishft(i,-1)
 119          !write(0,*) 'iand 2',i,ilog2_8
 120       endif
 121       end function ilog2_8
 122
 123       ! ----------------------------------------------------------------
 124
 125       function ilog2_4(i_in)
 126       implicit none
 127       integer(kind=4), value :: i_in
 128       integer(kind=4)::ilog2_4,i
 129       ilog2_4=0
 130       i=i_in
 131       if(i<=0) return
 132       if(iand(i,i-1)/=0) then
 133          !write(0,*) 'iand i-1'
 134          ilog2_4=1
 135       endif
 136       if(iand(i,Z'FFFF0000')/=0) then
 137          ilog2_4=ilog2_4+16
 138          i=ishft(i,-16)
 139          !write(0,*) 'iand ffff' ,i,ilog2_4
 140       endif
 141       if(iand(i,Z'0000FF00')/=0) then
 142          ilog2_4=ilog2_4+8
 143          i=ishft(i,-8)
 144          !write(0,*) 'iand ff',i,ilog2_4
 145       endif
 146       if(iand(i,Z'000000F0')/=0) then
 147          ilog2_4=ilog2_4+4
 148          i=ishft(i,-4)
 149          !write(0,*) 'iand f',i,ilog2_4
 150       endif
 151       if(iand(i,Z'0000000C')/=0) then
 152          ilog2_4=ilog2_4+2
 153          i=ishft(i,-2)
 154          !write(0,*) 'iand c',i,ilog2_4
 155       endif
 156       if(iand(i,Z'00000002')/=0) then
 157          ilog2_4=ilog2_4+1
 158          i=ishft(i,-1)
 159          !write(0,*) 'iand 2',i,ilog2_4
 160       endif
 161       end function ilog2_4
 162
 163       ! ----------------------------------------------------------------
 164
 165       function ilog2_2(i_in)
 166       implicit none
 167       integer(kind=2), value :: i_in
 168       integer(kind=2)::ilog2_2,i
 169       ilog2_2=0
 170       i=i_in
 171       if(i<=0) return
 172 ! WPS modification for the XL compiler
 173 !      if(iand(i,i-1)/=0) then
 174       if(iand(i,i-1_2)/=0) then
 175          !write(0,*) 'iand i-1'
 176          ilog2_2=1
 177       endif
 178       if(iand(i,Z'FF00')/=0) then
 179          ilog2_2=ilog2_2+8
 180          i=ishft(i,-8)
 181          !write(0,*) 'iand ff',i,ilog2_2
 182       endif
 183       if(iand(i,Z'00F0')/=0) then
 184          ilog2_2=ilog2_2+4
 185          i=ishft(i,-4)
 186          !write(0,*) 'iand f',i,ilog2_2
 187       endif
 188       if(iand(i,Z'000C')/=0) then
 189          ilog2_2=ilog2_2+2
 190          i=ishft(i,-2)
 191          !write(0,*) 'iand c',i,ilog2_2
 192       endif
 193       if(iand(i,Z'0002')/=0) then
 194          ilog2_2=ilog2_2+1
 195          i=ishft(i,-1)
 196          !write(0,*) 'iand 2',i,ilog2_2
 197       endif
 198       end function ilog2_2
 199
 200       ! ----------------------------------------------------------------
 201
 202       function ilog2_1(i_in)
 203       implicit none
 204       integer(kind=1), value :: i_in
 205       integer(kind=1)::ilog2_1,i
 206       ilog2_1=0
 207       i=i_in
 208       if(i<=0) return
 209 ! WPS modification for the XL compiler
 210 !      if(iand(i,i-1)/=0) then
 211       if(iand(i,i-1_1)/=0) then
 212          !write(0,*) 'iand i-1'
 213          ilog2_1=1
 214       endif
 215       if(iand(i,Z'F0')/=0) then
 216          ilog2_1=ilog2_1+4
 217          i=ishft(i,-4)
 218          !write(0,*) 'iand f',i,ilog2_1
 219       endif
 220       if(iand(i,Z'0C')/=0) then
 221          ilog2_1=ilog2_1+2
 222          i=ishft(i,-2)
 223          !write(0,*) 'iand c',i,ilog2_1
 224       endif
 225       if(iand(i,Z'02')/=0) then
 226          ilog2_1=ilog2_1+1
 227          i=ishft(i,-1)
 228          !write(0,*) 'iand 2',i,ilog2_1
 229       endif
 230       end function ilog2_1
 231       end module intmath
 232 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 233 c$$$ TEST PROGRAM FOR THIS MODULE
 234 c$$$      program test_intmath
 235 c$$$      use intmath
 236 c$$$      implicit none
 237 c$$$      real(kind=16) :: temp
 238 c$$$      real(kind=16), parameter :: alog2=log(2.0_16)
 239 c$$$      integer(kind=8), parameter  ::                                    &
 240 c$$$     &     one=1,big=Z'7FFFFFFFFFFFFFFF',small=-2000000_8,              &
 241 c$$$     &     check=Z'1FFFFFFF'
 242 c$$$      integer(kind=8) :: ival, iret
 243 c$$$      !$OMP PARALLEL DO PRIVATE(ival,temp,iret)
 244 c$$$      do ival=small,big
 245 c$$$ 10      format(Z16,' -- MISMATCH: ',I0,'=>',I0,' (',I0,' = ',F0.10,')')
 246 c$$$ 20      format(Z16,' -- OKAY:     ',I0,'=>',I0,' (',I0,' = ',F0.10,')')
 247 c$$$         if(ival+one<ival) then
 248 c$$$            temp=log(real(max(ival,one),kind=16))/alog2
 249 c$$$         else
 250 c$$$            temp=log(real(max(ival+one,one),kind=16))/alog2
 251 c$$$         endif
 252 c$$$         iret=i1log2(ival)
 253 c$$$         if(iret/=ceiling(temp) .or. ival==0 .or. ival==check) then
 254 c$$$            !$OMP CRITICAL
 255 c$$$            if(iret/=ceiling(temp)) then
 256 c$$$               print 10, ival, ival, iret,ceiling(temp),temp
 257 c$$$            else
 258 c$$$               print 20, ival, ival, iret,ceiling(temp),temp
 259 c$$$            endif
 260 c$$$            !$OMP END CRITICAL
 261 c$$$         endif
 262 c$$$      enddo
 263 c$$$      !$OMP END PARALLEL DO
 264 c$$$      end program test_intmath