!python3 -m numpy.f2py -c -m BIOVAR biovar.f90

!19 Bioclimate variables defined in USGS ds691.pdf
!https://pubs.usgs.gov/ds/691/
subroutine biovar(t2m,tmax,tmin,rain,ix,iy,it,&
                  b1,b2,b3,b4,b4a,b5,b6,b7,b8,&
                  b9,b10,b11,b12,b13,b14,b15, &
                  b16,b17,b18,b19)
!F2PY INTENT(IN)::t2m,tmax,tmin,rain
!F2PY INTENT(OUT)::b1,b2,b3,b4,b4a,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19
!F2PY INTENT(HIDE)::ix,iy,it
integer ix           !number of longitudes
integer iy           !number of latitudes
integer it           !number of months, 14 starting from Jan
real t2m(it,iy,ix)   !monthly mean temperature, Celcius
real tmax(it,iy,ix)  !monthly maximum temperature, Celcius
real tmin(it,iy,ix)  !monthly minimum temperature, Celcius
real rain(it,iy,ix)  !monthly total precipitation, mm
real b1(iy,ix),b2(iy,ix),b3(iy,ix),b4(iy,ix),b4a(iy,ix),b5(iy,ix)     
real b6(iy,ix),b7(iy,ix),b8(iy,ix),b9(iy,ix),b10(iy,ix),b11(iy,ix)
real b12(iy,ix),b13(iy,ix),b14(iy,ix),b15(iy,ix),b16(iy,ix)     
real b17(iy,ix),b18(iy,ix),b19(iy,ix)     
! local variables
integer i,j,k
integer it2
real s1,s2,s3,s4,tmp
integer k1,k2,k3,k4

it2=it-2
do j=1,iy
do i=1,ix
 b1(j,i)=0.
 b2(j,i)=0.
 b5(j,i)=-1.e10
 b6(j,i)=1.e10
 b12(j,i)=0.
 b13(j,i)=-1.e10
 b14(j,i)=1.e10
 do k=1,it2
  b1(j,i)=b1(j,i)+t2m(k,j,i)
  b2(j,i)=b2(j,i)+tmax(k,j,i)-tmin(k,j,i)
  if(b5(j,i) .le. tmax(k,j,i))b5(j,i)=tmax(k,j,i)
  if(b6(j,i) .ge. tmin(k,j,i))b6(j,i)=tmin(k,j,i)
  b12(j,i)=b12(j,i)+rain(k,j,i)
  if(b13(j,i) .le. rain(k,j,i))b13(j,i)=rain(k,j,i)
  if(b14(j,i) .ge. rain(k,j,i))b14(j,i)=rain(k,j,i)
 enddo
 b1(j,i)=b1(j,i)/it2
 b2(j,i)=b2(j,i)/it2
 b7(j,i)=b5(j,i)-b6(j,i)
 b3(j,i)=b2(j,i)/b7(j,i)*100
enddo
enddo

do j=1,iy
do i=1,ix
 b4(j,i)=0.
 b15(j,i)=0.
 do k=1,it2
  b4(j,i)=b4(j,i)+(t2m(k,j,i)-b1(j,i))**2
  b15(j,i)=b15(j,i)+(rain(k,j,i)-b12(j,i)/it2)**2
 enddo
 b4(j,i)=sqrt(b4(j,i)/(it2-1))
 b15(j,i)=sqrt(b15(j,i)/(it2-1))
 b15(j,i)=b15(j,i)/(1+b12(j,i)/12)*100
 b4a(j,i)=b4(j,i)/(b1(j,i)+273.15)*100
enddo
enddo

do j=1,iy
do i=1,ix
 s1=-1.e10
 s2=1.e10
 s3=-1.e10
 s4=1.e10
 do k=1,it2
  tmp=rain(k,j,i)+rain(k+1,j,i)+rain(k+2,j,i)
  if(s1 .le. tmp) then
   s1=tmp
   b16(j,i)=s1
   k1=k
  endif
  if(s2 .ge. tmp) then
   k2=k
   s2=tmp
   b17(j,i)=s2
  endif
  tmp=t2m(k,j,i)+t2m(k+1,j,i)+t2m(k+2,j,i)
  if(s3 .le. tmp) then
   k3=k
   s3=tmp
  endif
  if(s4 .ge. tmp) then
   k4=k
   s4=tmp
  endif
 enddo
 b8(j,i)=(t2m(k1,j,i)+t2m(k1+1,j,i)+t2m(k1+2,j,i))/3
 b18(j,i)=(rain(k1,j,i)+rain(k1+1,j,i)+rain(k1+2,j,i))/3
 b9(j,i)=(t2m(k2,j,i)+t2m(k2+1,j,i)+t2m(k2+2,j,i))/3
 b19(j,i)=(rain(k2,j,i)+rain(k2+1,j,i)+rain(k2+2,j,i))/3
 b10(j,i)=(t2m(k3,j,i)+t2m(k3+1,j,i)+t2m(k3+2,j,i))/3
 b11(j,i)=(t2m(k4,j,i)+t2m(k4+1,j,i)+t2m(k4+2,j,i))/3
enddo
enddo

return
end subroutine biovar