module eq_regions_mod

!

!     Purpose.

!     --------

!           eq_regions_mod provides the code to perform a high level

!           partitioning of the surface of a sphere into regions of

!           equal area and small diameter.

!           the type.

!

!     Background.

!     -----------

!     This Fortran version of eq_regions is a much cut down version of the

!     "Recursive Zonal Equal Area (EQ) Sphere Partitioning Toolbox" of the

!     same name developed by Paul Leopardi at the University of New South Wales.

!     This version has been coded specifically for the case of partitioning the

!     surface of a sphere or S^dim (where dim=2) as denoted in the original code.

!     Only a subset of the original eq_regions package has been coded to determine

!     the high level distribution of regions on a sphere, as the detailed

!     distribution of grid points to each region is left to IFS software.

!     This is required to take into account the spatial distribution of grid

!     points in an IFS gaussian grid and provide an optimal (i.e. exact)

!     distribution of grid points over regions.

!

!     The following copyright notice for the eq_regions package is included from

!     the original MatLab release.

!

!     +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!     + Release 1.10 2005-06-26                                                 +

!     +                                                                         +

!     + Copyright (c) 2004, 2005, University of New South Wales                 +

!     +                                                                         +

!     + Permission is hereby granted, free of charge, to any person obtaining   +

!     + a copy of this software and associated documentation files (the         +

!     + "Software"), to deal in the Software without restriction, including     +

!     + without limitation the rights to use, copy, modify, merge, publish,     +

!     + distribute, sublicense, and/or sell copies of the Software, and to      +

!     + permit persons to whom the Software is furnished to do so, subject to   +

!     + the following conditions:                                               +

!     +                                                                         +

!     + The above copyright notice and this permission notice shall be included +

!     + in all copies or substantial portions of the Software.                  +

!     +                                                                         +

!     + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         +

!     + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      +

!     + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  +

!     + IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY    +

!     + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,    +

!     + TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE       +

!     + SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                  +

!     +                                                                         +

!     +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!

!     Author.

!     -------

!        George Mozdzynski *ECMWF*

!

!     Modifications.

!     --------------

!        Original : 2006-04-15

!

!--------------------------------------------------------------------------------

USE PARKIND1  ,ONLY : JPIM     ,JPRB

IMPLICIT NONE

SAVE

PRIVATE

PUBLIC eq_regions,l_regions_debug,n_regions_ns,n_regions_ew,n_regions,my_region_ns,my_region_ew

real(kind=jprb) pi

logical :: l_regions_debug=.false.

integer(kind=jpim) :: n_regions_ns

integer(kind=jpim) :: n_regions_ew

integer(kind=jpim) :: my_region_ns

integer(kind=jpim) :: my_region_ew

integer(kind=jpim),allocatable :: n_regions(:)

!$OMP THREADPRIVATE(l_regions_debug,my_region_ew,my_region_ns,n_regions_ew,n_regions_ns,pi,n_regions)

CONTAINS

!

! eq_regions uses the zonal equal area sphere partitioning algorithm to partition

! the surface of a sphere into N regions of equal area and small diameter.

!

integer(kind=jpim),intent(in) :: N

integer(kind=jpim) :: n_collars,j

real(kind=jprb),allocatable :: r_regions(:)

real(kind=jprb) :: c_polar

  !

  ! We have only one region, which must be the whole sphere.

  !

else

  !

  ! Given N, determine c_polar

  ! the colatitude of the North polar spherical cap.

  !

  !

  ! Given N, determine the ideal angle for spherical collars.

  ! Based on N, this ideal angle, and c_polar,

  ! determine n_collars, the number of collars between the polar caps.

  !

  !

  ! Given N, c_polar and n_collars, determine r_regions,

  ! a list of the ideal real number of regions in each collar,

  ! plus the polar caps.

  ! The number of elements is n_collars+2.

  ! r_regions[1] is 1.

  ! r_regions[n_collars+2] is 1.

  ! The sum of r_regions is N.

  !

  ! Given N and r_regions, determine n_regions, a list of the natural number

  ! of regions in each collar and the polar caps.

  ! This list is as close as possible to r_regions.

  ! The number of elements is n_collars+2.

  ! n_regions[1] is 1.

  ! n_regions[n_collars+2] is 1.

  ! The sum of n_regions is N.

  !

  endif

endif

  enddo

endif

end subroutine eq_regions

!

!NUM_COLLARS The number of collars between the polar caps

!

! Given N, an ideal angle, and c_polar,

! determine n_collars, the number of collars between the polar caps.

!

integer(kind=jpim),intent(in) :: N

real(kind=jprb),intent(in) :: a_ideal,c_polar

integer(kind=jpim) :: num_c

logical enough

if( enough )then

else

  num_c = 0

endif

return

end function num_collars

!

!IDEAL_REGION_LIST The ideal real number of regions in each zone

!

! List the ideal real number of regions in each collar, plus the polar caps.

!

! Given N, c_polar and n_collars, determine r_regions, a list of the ideal real

! number of regions in each collar, plus the polar caps.

! The number of elements is n_collars+2.

! r_regions[1] is 1.

! r_regions[n_collars+2] is 1.

! The sum of r_regions is N.

!

integer(kind=jpim),intent(in) :: N,n_collars

real(kind=jprb),intent(in) :: c_polar

real(kind=jprb),intent(out) :: r_regions(n_collars+2)

integer(kind=jpim) :: collar_n

real(kind=jprb) :: ideal_region_area,ideal_collar_area

real(kind=jprb) :: a_fitting

  !

  ! Based on n_collars and c_polar, determine a_fitting,

  ! the collar angle such that n_collars collars fit between the polar caps.

  !

    ideal_collar_area = area_of_collar(c_polar+(collar_n-1)*a_fitting, &

  enddo

endif

end subroutine ideal_region_list

function ideal_collar_angle(N) result(ideal)

!

! IDEAL_COLLAR_ANGLE The ideal angle for spherical collars of an EQ partition

!

! IDEAL_COLLAR_ANGLE(N) sets ANGLE to the ideal angle for the

! spherical collars of an EQ partition of the unit sphere S^2 into N regions.

!

integer(kind=jpim),intent(in) :: N

real(kind=jprb) :: ideal

return

end function ideal_collar_angle

!

! ROUND_TO_NATURALS Round off a given list of numbers of regions

!

! Given N and r_regions, determine n_regions, a list of the natural number

! of regions in each collar and the polar caps.

! This list is as close as possible to r_regions, using rounding.

! The number of elements is n_collars+2.

! n_regions[1] is 1.

! n_regions[n_collars+2] is 1.

! The sum of n_regions is N.

!

integer(kind=jpim),intent(in) :: N,n_collars

real(kind=jprb),intent(in) :: r_regions(n_collars+2)

integer(kind=jpim) :: zone_n

real(kind=jprb) :: discrepancy

discrepancy = 0.0_jprb

enddo

end subroutine round_to_naturals

!

! Given N, determine the colatitude of the North polar spherical cap.

!

integer(kind=jpim),intent(in) :: N

real(kind=jprb) :: area

real(kind=jprb) :: polar_c

endif

return

end function polar_colat

!

! AREA_OF_IDEAL_REGION(N) sets AREA to be the area of one of N equal

! area regions on S^2, that is 1/N times AREA_OF_SPHERE.

!

integer(kind=jpim),intent(in) :: N

real(kind=jprb) :: area_of_sphere

real(kind=jprb) :: area

return

end function area_of_ideal_region

!

! SRADIUS_OF_CAP(AREA) returns the spherical radius of

! an S^2 spherical cap of area AREA.

!

real(kind=jprb),intent(in) :: area

real(kind=jprb) :: sradius

return

end function sradius_of_cap

!

! AREA_OF_COLLAR Area of spherical collar

!

! AREA_OF_COLLAR(A_TOP, A_BOT) sets AREA to be the area of an S^2 spherical

! collar specified by A_TOP, A_BOT, where A_TOP is top (smaller) spherical radius,

! A_BOT is bottom (larger) spherical radius.

!

real(kind=jprb),intent(in) :: a_top,a_bot

real(kind=jprb) area

return

end function area_of_collar

function area_of_cap(s_cap) result(area)

!

! AREA_OF_CAP Area of spherical cap

!

! AREA_OF_CAP(S_CAP) sets AREA to be the area of an S^2 spherical

! cap of spherical radius S_CAP.

!

real(kind=jprb),intent(in) :: s_cap

real(kind=jprb) area

return

end function area_of_cap

real(kind=jprb),intent(in) :: x

real(kind=jprb) :: gamma_res

real(kind=jprb) :: p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13

real(kind=jprb) :: w,y

integer(kind=jpim) :: k,n

parameter (&

& p0 =   0.999999999999999990e+00_jprb,&

& p1 =  -0.422784335098466784e+00_jprb,&

& p2 =  -0.233093736421782878e+00_jprb,&

& p3 =   0.191091101387638410e+00_jprb,&

& p4 =  -0.024552490005641278e+00_jprb,&

& p5 =  -0.017645244547851414e+00_jprb,&

& p6 =   0.008023273027855346e+00_jprb)

parameter (&

& p7 =  -0.000804329819255744e+00_jprb,&

& p8 =  -0.000360837876648255e+00_jprb,&

& p9 =   0.000145596568617526e+00_jprb,&

& p10 = -0.000017545539395205e+00_jprb,&

& p11 = -0.000002591225267689e+00_jprb,&

& p12 =  0.000001337767384067e+00_jprb,&

& p13 = -0.000000199542863674e+00_jprb)

y = ((((((((((((p13 * w + p12) * w + p11) * w + p10) *&

&    w + p9) * w + p8) * w + p7) * w + p6) * w + p5) *&

  end do

else

  w = 1

  end do

end if

return

end function gamma

end module eq_regions_mod