diff --git a/.fprettify.rc b/.fprettify.rc new file mode 100644 index 0000000..27861c9 --- /dev/null +++ b/.fprettify.rc @@ -0,0 +1,19 @@ +enable-replacements=True +c-relations=True +case=[1,1,1,1] +indent=4 +line-length=132 +strict-indent=True +strip-comments=True +whitespace-relational=True +whitespace-logical=True +whitespace-plusminus=True +whitespace-multdiv=True +whitespace-comma=True +whitespace-intrinsics=True +whitespace-print=False +whitespace-type=False +whitespace-decl=True +enable-decl=True +disable-fypp=True +disable-indent-mod=False diff --git a/radbelt.code-workspace b/radbelt.code-workspace index 489f23c..7f2c28f 100644 --- a/radbelt.code-workspace +++ b/radbelt.code-workspace @@ -8,6 +8,11 @@ "files.trimTrailingWhitespace": true, "editor.insertSpaces": true, "editor.tabSize": 4, - "editor.trimAutoWhitespace": true + "editor.trimAutoWhitespace": true, + "fortran.linter.includePaths": [ + "${workspaceFolder}/src", + "${workspaceFolder}/test" + ], + "fortran.linter.initialize": false } } \ No newline at end of file diff --git a/src/radbelt_c_module.f90 b/src/radbelt_c_module.f90 index 94c7a21..e81fc1a 100644 --- a/src/radbelt_c_module.f90 +++ b/src/radbelt_c_module.f90 @@ -2,7 +2,7 @@ !> ! Experimental C interface to the radbelt module. - module radbelt_c_module +module radbelt_c_module use iso_c_binding, only: c_double, c_int, c_char, c_null_char, & c_intptr_t, c_ptr, c_loc, c_f_pointer, & @@ -11,194 +11,194 @@ module radbelt_c_module implicit none - contains +contains !***************************************************************************************** !***************************************************************************************** !> ! Convert C string to Fortran -function c2f_str(cstr) result(fstr) + function c2f_str(cstr) result(fstr) - character(kind=c_char,len=1),dimension(:),intent(in) :: cstr !! string from C - character(len=:),allocatable :: fstr !! fortran string + character(kind=c_char, len=1), dimension(:), intent(in) :: cstr !! string from C + character(len=:), allocatable :: fstr !! fortran string - integer :: i !! counter + integer :: i !! counter - fstr = '' - do i = 1, size(cstr) - fstr = fstr//cstr(i) - end do - fstr = trim(fstr) + fstr = '' + do i = 1, size(cstr) + fstr = fstr//cstr(i) + end do + fstr = trim(fstr) -end function c2f_str + end function c2f_str !***************************************************************************************** !> ! Convert an integer pointer to a [[radbelt_type]] pointer. -subroutine int_pointer_to_f_pointer(ipointer, p) + subroutine int_pointer_to_f_pointer(ipointer, p) - integer(c_intptr_t),intent(in) :: ipointer !! integer pointer from C - type(radbelt_type),pointer :: p !! fortran pointer + integer(c_intptr_t), intent(in) :: ipointer !! integer pointer from C + type(radbelt_type), pointer :: p !! fortran pointer - type(c_ptr) :: cp + type(c_ptr) :: cp - cp = transfer(ipointer, c_null_ptr) - if (c_associated(cp)) then - call c_f_pointer(cp, p) - else - p => null() - end if + cp = transfer(ipointer, c_null_ptr) + if (c_associated(cp)) then + call c_f_pointer(cp, p) + else + p => null() + end if -end subroutine int_pointer_to_f_pointer + end subroutine int_pointer_to_f_pointer !***************************************************************************************** !> ! create a [[radbelt_type]] from C -subroutine initialize_c(ipointer) bind(C, name="initialize_c") + subroutine initialize_c(ipointer) bind(C, name="initialize_c") - integer(c_intptr_t),intent(out) :: ipointer - type(radbelt_type),pointer :: p - type(c_ptr) :: cp + integer(c_intptr_t), intent(out) :: ipointer + type(radbelt_type), pointer :: p + type(c_ptr) :: cp - allocate(p) - cp = c_loc(p) - ipointer = transfer(cp, 0_c_intptr_t) + allocate (p) + cp = c_loc(p) + ipointer = transfer(cp, 0_c_intptr_t) -end subroutine initialize_c + end subroutine initialize_c !***************************************************************************************** !> ! destroy a [[radbelt_type]] from C -subroutine destroy_c(ipointer) bind(C, name="destroy_c") + subroutine destroy_c(ipointer) bind(C, name="destroy_c") - integer(c_intptr_t),intent(in) :: ipointer - type(radbelt_type),pointer :: p + integer(c_intptr_t), intent(in) :: ipointer + type(radbelt_type), pointer :: p - call int_pointer_to_f_pointer(ipointer,p) - if (associated(p)) deallocate(p) + call int_pointer_to_f_pointer(ipointer, p) + if (associated(p)) deallocate (p) -end subroutine destroy_c + end subroutine destroy_c !***************************************************************************************** !> ! C interface for setting the `trm` data file path -subroutine set_trm_file_path_c(ipointer, aep8_dir, n) bind(C, name="set_trm_file_path_c") + subroutine set_trm_file_path_c(ipointer, aep8_dir, n) bind(C, name="set_trm_file_path_c") - integer(c_intptr_t),intent(in) :: ipointer - integer(c_int),intent(in) :: n !! size of `aep8_dir` - character(kind=c_char,len=1),dimension(n),intent(in) :: aep8_dir + integer(c_intptr_t), intent(in) :: ipointer + integer(c_int), intent(in) :: n !! size of `aep8_dir` + character(kind=c_char, len=1), dimension(n), intent(in) :: aep8_dir - character(len=:),allocatable :: aep8_dir_ - type(radbelt_type),pointer :: p + character(len=:), allocatable :: aep8_dir_ + type(radbelt_type), pointer :: p - call int_pointer_to_f_pointer(ipointer, p) + call int_pointer_to_f_pointer(ipointer, p) - if (associated(p)) then - aep8_dir_ = c2f_str(aep8_dir) - call p%set_trm_file_path(aep8_dir_) - else - error stop 'error in set_trm_file_path_c: class is not associated' - end if + if (associated(p)) then + aep8_dir_ = c2f_str(aep8_dir) + call p%set_trm_file_path(aep8_dir_) + else + error stop 'error in set_trm_file_path_c: class is not associated' + end if - end subroutine set_trm_file_path_c + end subroutine set_trm_file_path_c !***************************************************************************************** !***************************************************************************************** !> ! C interface for setting the `igrf` data file path - subroutine set_igrf_file_path_c(ipointer, igrf_dir, n) bind(C, name="set_igrf_file_path") + subroutine set_igrf_file_path_c(ipointer, igrf_dir, n) bind(C, name="set_igrf_file_path") - integer(c_intptr_t),intent(in) :: ipointer - integer(c_int),intent(in) :: n !! size of `igrf_dir` - character(kind=c_char,len=1),dimension(n),intent(in) :: igrf_dir + integer(c_intptr_t), intent(in) :: ipointer + integer(c_int), intent(in) :: n !! size of `igrf_dir` + character(kind=c_char, len=1), dimension(n), intent(in) :: igrf_dir - character(len=:),allocatable :: igrf_dir_ - type(radbelt_type),pointer :: p + character(len=:), allocatable :: igrf_dir_ + type(radbelt_type), pointer :: p - call int_pointer_to_f_pointer(ipointer, p) + call int_pointer_to_f_pointer(ipointer, p) - if (associated(p)) then - igrf_dir_ = c2f_str(igrf_dir) - call p%set_igrf_file_path(igrf_dir_) - else - error stop 'error in set_igrf_file_path: class is not associated' - end if + if (associated(p)) then + igrf_dir_ = c2f_str(igrf_dir) + call p%set_igrf_file_path(igrf_dir_) + else + error stop 'error in set_igrf_file_path: class is not associated' + end if - end subroutine set_igrf_file_path_c + end subroutine set_igrf_file_path_c !***************************************************************************************** !***************************************************************************************** !> ! C interface for setting the data file paths - subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir, n, m) bind(C, name="set_data_files_paths_c") + subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir, n, m) bind(C, name="set_data_files_paths_c") - integer(c_intptr_t),intent(in) :: ipointer - integer(c_int),intent(in) :: n !! size of `aep8_dir` - character(kind=c_char,len=1),dimension(n),intent(in) :: aep8_dir - integer(c_int),intent(in) :: m !! size of `igrf_dir` - character(kind=c_char,len=1),dimension(m),intent(in) :: igrf_dir + integer(c_intptr_t), intent(in) :: ipointer + integer(c_int), intent(in) :: n !! size of `aep8_dir` + character(kind=c_char, len=1), dimension(n), intent(in) :: aep8_dir + integer(c_int), intent(in) :: m !! size of `igrf_dir` + character(kind=c_char, len=1), dimension(m), intent(in) :: igrf_dir - character(len=:),allocatable :: aep8_dir_, igrf_dir_ - type(radbelt_type),pointer :: p + character(len=:), allocatable :: aep8_dir_, igrf_dir_ + type(radbelt_type), pointer :: p - call int_pointer_to_f_pointer(ipointer, p) + call int_pointer_to_f_pointer(ipointer, p) - if (associated(p)) then + if (associated(p)) then - aep8_dir_ = c2f_str(aep8_dir) - igrf_dir_ = c2f_str(igrf_dir) + aep8_dir_ = c2f_str(aep8_dir) + igrf_dir_ = c2f_str(igrf_dir) - call p%set_data_files_paths(aep8_dir_, igrf_dir_) + call p%set_data_files_paths(aep8_dir_, igrf_dir_) - else - error stop 'error in set_data_files_paths_c: class is not associated' - end if + else + error stop 'error in set_data_files_paths_c: class is not associated' + end if - end subroutine set_data_files_paths_c + end subroutine set_data_files_paths_c !***************************************************************************************** !***************************************************************************************** !> ! C interface to [[get_flux_g]]. -subroutine get_flux_g_c(ipointer,lon,lat,height,year,e,imname,flux) bind(C, name="get_flux_g_c") + subroutine get_flux_g_c(ipointer, lon, lat, height, year, e, imname, flux) bind(C, name="get_flux_g_c") - integer(c_intptr_t),intent(in) :: ipointer - real(c_double),intent(in) :: lon !! geodetic longitude in degrees (east) - real(c_double),intent(in) :: lat !! geodetic latitude in degrees (north) - real(c_double),intent(in) :: height !! altitude in km above sea level - real(c_double),intent(in) :: year !! decimal year for which geomagnetic field is to + integer(c_intptr_t), intent(in) :: ipointer + real(c_double), intent(in) :: lon !! geodetic longitude in degrees (east) + real(c_double), intent(in) :: lat !! geodetic latitude in degrees (north) + real(c_double), intent(in) :: height !! altitude in km above sea level + real(c_double), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) - real(c_double),intent(in) :: e !! minimum energy - integer(c_int),intent(in) :: imname !! which method to use: + real(c_double), intent(in) :: e !! minimum energy + integer(c_int), intent(in) :: imname !! which method to use: !! !! * 1 -- particle species: electrons, solar activity: min !! * 2 -- particle species: electrons, solar activity: max !! * 3 -- particle species: protons, solar activity: min !! * 4 -- particle species: protons, solar activity: max - real(c_double),intent(out) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. + real(c_double), intent(out) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. - type(radbelt_type),pointer :: p + type(radbelt_type), pointer :: p - call int_pointer_to_f_pointer(ipointer, p) + call int_pointer_to_f_pointer(ipointer, p) - if (associated(p)) then + if (associated(p)) then - flux = p%get_flux(lon,lat,height,year,e,imname) + flux = p%get_flux(lon, lat, height, year, e, imname) - else - error stop 'error in get_flux_g_c: class is not associated' - end if + else + error stop 'error in get_flux_g_c: class is not associated' + end if -end subroutine get_flux_g_c + end subroutine get_flux_g_c !***************************************************************************************** - end module radbelt_c_module +end module radbelt_c_module !***************************************************************************************** \ No newline at end of file diff --git a/src/radbelt_kinds_module.F90 b/src/radbelt_kinds_module.F90 index b814007..de741b7 100644 --- a/src/radbelt_kinds_module.F90 +++ b/src/radbelt_kinds_module.F90 @@ -2,36 +2,36 @@ !> ! Numeric kind definitions for radbelt. - module radbelt_kinds_module +module radbelt_kinds_module - use,intrinsic :: iso_fortran_env + use, intrinsic :: iso_fortran_env implicit none private #ifdef REAL32 - integer,parameter,public :: wp = real32 !! Real working precision [4 bytes] + integer, parameter, public :: wp = real32 !! Real working precision [4 bytes] #elif REAL64 - integer,parameter,public :: wp = real64 !! Real working precision [8 bytes] + integer, parameter, public :: wp = real64 !! Real working precision [8 bytes] #elif REAL128 - integer,parameter,public :: wp = real128 !! Real working precision [16 bytes] + integer, parameter, public :: wp = real128 !! Real working precision [16 bytes] #else - integer,parameter,public :: wp = real64 !! Real working precision if not specified [8 bytes] + integer, parameter, public :: wp = real64 !! Real working precision if not specified [8 bytes] #endif #ifdef INT8 - integer,parameter,public :: ip = int8 !! Integer working precision [1 byte] + integer, parameter, public :: ip = int8 !! Integer working precision [1 byte] #elif INT16 - integer,parameter,public :: ip = int16 !! Integer working precision [2 bytes] + integer, parameter, public :: ip = int16 !! Integer working precision [2 bytes] #elif INT32 - integer,parameter,public :: ip = int32 !! Integer working precision [4 bytes] + integer, parameter, public :: ip = int32 !! Integer working precision [4 bytes] #elif INT64 - integer,parameter,public :: ip = int64 !! Integer working precision [8 bytes] + integer, parameter, public :: ip = int64 !! Integer working precision [8 bytes] #else - integer,parameter,public :: ip = int32 !! Integer working precision if not specified [4 bytes] + integer, parameter, public :: ip = int32 !! Integer working precision if not specified [4 bytes] #endif !***************************************************************************************** - end module radbelt_kinds_module +end module radbelt_kinds_module !***************************************************************************************** diff --git a/src/radbelt_module.f90 b/src/radbelt_module.f90 index 7a0ff78..2f30265 100644 --- a/src/radbelt_module.f90 +++ b/src/radbelt_module.f90 @@ -8,61 +8,61 @@ module radbelt_module - use radbelt_kinds_module - use trmfun_module - use shellig_module - - implicit none - - type,public :: radbelt_type - !! the main class that can be used to get the flux. - private - type(trm_type) :: trm - type(shellig_type) :: igrf - contains - private - generic,public :: get_flux => get_flux_g_, get_flux_c_ - procedure :: get_flux_g_, get_flux_c_ - procedure,public :: set_trm_file_path, & - set_igrf_file_path, & - set_data_files_paths - end type radbelt_type - - interface get_flux - !! simple function versions for testing - procedure :: get_flux_g - procedure :: get_flux_c - end interface - public :: get_flux - - contains + use radbelt_kinds_module + use trmfun_module + use shellig_module + + implicit none + + type, public :: radbelt_type + !! the main class that can be used to get the flux. + private + type(trm_type) :: trm + type(shellig_type) :: igrf + contains + private + generic, public :: get_flux => get_flux_g_, get_flux_c_ + procedure :: get_flux_g_, get_flux_c_ + procedure, public :: set_trm_file_path, & + set_igrf_file_path, & + set_data_files_paths + end type radbelt_type + + interface get_flux + !! simple function versions for testing + procedure :: get_flux_g + procedure :: get_flux_c + end interface + public :: get_flux + +contains !***************************************************************************************** !> ! Set the `trm` path. - subroutine set_trm_file_path(me, dir) + subroutine set_trm_file_path(me, dir) - class(radbelt_type),intent(inout) :: me - character(len=*),intent(in) :: dir + class(radbelt_type), intent(inout) :: me + character(len=*), intent(in) :: dir - call me%trm%set_data_file_dir(trim(dir)) + call me%trm%set_data_file_dir(trim(dir)) - end subroutine set_trm_file_path + end subroutine set_trm_file_path !***************************************************************************************** !***************************************************************************************** !> ! Set the `igrf` path. - subroutine set_igrf_file_path(me, dir) + subroutine set_igrf_file_path(me, dir) - class(radbelt_type),intent(inout) :: me - character(len=*),intent(in) :: dir + class(radbelt_type), intent(inout) :: me + character(len=*), intent(in) :: dir - call me%igrf%set_data_file_dir(trim(dir)) + call me%igrf%set_data_file_dir(trim(dir)) - end subroutine set_igrf_file_path + end subroutine set_igrf_file_path !***************************************************************************************** !***************************************************************************************** @@ -71,46 +71,46 @@ end subroutine set_igrf_file_path ! If not used or blank, the folder `data/aep8` and `data/igrf` in the ! current working directory is assumed - subroutine set_data_files_paths(me, aep8_dir, igrf_dir) + subroutine set_data_files_paths(me, aep8_dir, igrf_dir) - class(radbelt_type),intent(inout) :: me - character(len=*),intent(in) :: aep8_dir - character(len=*),intent(in) :: igrf_dir + class(radbelt_type), intent(inout) :: me + character(len=*), intent(in) :: aep8_dir + character(len=*), intent(in) :: igrf_dir - call me%set_trm_file_path(trim(aep8_dir)) - call me%set_igrf_file_path(trim(igrf_dir)) + call me%set_trm_file_path(trim(aep8_dir)) + call me%set_igrf_file_path(trim(igrf_dir)) - end subroutine set_data_files_paths + end subroutine set_data_files_paths !***************************************************************************************** !***************************************************************************************** !> ! Calculate the flux of trapped particles at a specific location and time. - function get_flux_g_(me,lon,lat,height,year,e,imname) result(flux) + function get_flux_g_(me, lon, lat, height, year, e, imname) result(flux) - class(radbelt_type),intent(inout) :: me - real(wp),intent(in) :: lon !! geodetic longitude in degrees (east) - real(wp),intent(in) :: lat !! geodetic latitude in degrees (north) - real(wp),intent(in) :: height !! altitude in km above sea level - real(wp),intent(in) :: year !! decimal year for which geomagnetic field is to + class(radbelt_type), intent(inout) :: me + real(wp), intent(in) :: lon !! geodetic longitude in degrees (east) + real(wp), intent(in) :: lat !! geodetic latitude in degrees (north) + real(wp), intent(in) :: height !! altitude in km above sea level + real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) - real(wp),intent(in) :: e !! minimum energy - integer,intent(in) :: imname !! which method to use: + real(wp), intent(in) :: e !! minimum energy + integer, intent(in) :: imname !! which method to use: !! !! * 1 -- particle species: electrons, solar activity: min !! * 2 -- particle species: electrons, solar activity: max !! * 3 -- particle species: protons, solar activity: min !! * 4 -- particle species: protons, solar activity: max - real(wp) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. + real(wp) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. - real(wp) :: xl !! l value - real(wp) :: bbx + real(wp) :: xl !! l value + real(wp) :: bbx - call me%igrf%igrf(lon,lat,height,year,xl,bbx) - call me%trm%aep8(e,xl,bbx,imname,flux) + call me%igrf%igrf(lon, lat, height, year, xl, bbx) + call me%trm%aep8(e, xl, bbx, imname, flux) - end function get_flux_g_ + end function get_flux_g_ !***************************************************************************************** !***************************************************************************************** @@ -121,27 +121,27 @@ end function get_flux_g_ !@note This routine is not efficient at all since it will reload all the ! files every time it is called. - function get_flux_g(lon,lat,height,year,e,imname) result(flux) + function get_flux_g(lon, lat, height, year, e, imname) result(flux) - real(wp),intent(in) :: lon !! geodetic longitude in degrees (east) - real(wp),intent(in) :: lat !! geodetic latitude in degrees (north) - real(wp),intent(in) :: height !! altitude in km above sea level - real(wp),intent(in) :: year !! decimal year for which geomagnetic field is to + real(wp), intent(in) :: lon !! geodetic longitude in degrees (east) + real(wp), intent(in) :: lat !! geodetic latitude in degrees (north) + real(wp), intent(in) :: height !! altitude in km above sea level + real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) - real(wp),intent(in) :: e !! minimum energy - integer,intent(in) :: imname !! which method to use: + real(wp), intent(in) :: e !! minimum energy + integer, intent(in) :: imname !! which method to use: !! !! * 1 -- particle species: electrons, solar activity: min !! * 2 -- particle species: electrons, solar activity: max !! * 3 -- particle species: protons, solar activity: min !! * 4 -- particle species: protons, solar activity: max - real(wp) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. + real(wp) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. - type(radbelt_type) :: radbelt + type(radbelt_type) :: radbelt - flux = radbelt%get_flux(lon,lat,height,year,e,imname) + flux = radbelt%get_flux(lon, lat, height, year, e, imname) - end function get_flux_g + end function get_flux_g !***************************************************************************************** !***************************************************************************************** @@ -149,28 +149,28 @@ end function get_flux_g ! Calculate the flux of trapped particles at a specific location and time. ! This is an alternate version of [[get_flux_g_]] for cartesian coordinates. - function get_flux_c_(me,v,year,e,imname) result(flux) + function get_flux_c_(me, v, year, e, imname) result(flux) - class(radbelt_type),intent(inout) :: me - real(wp),dimension(3),intent(in) :: v - real(wp),intent(in) :: year !! decimal year for which geomagnetic field is to + class(radbelt_type), intent(inout) :: me + real(wp), dimension(3), intent(in) :: v + real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) - real(wp),intent(in) :: e !! minimum energy - integer,intent(in) :: imname !! which method to use: + real(wp), intent(in) :: e !! minimum energy + integer, intent(in) :: imname !! which method to use: !! !! * 1 -- particle species: electrons, solar activity: min !! * 2 -- particle species: electrons, solar activity: max !! * 3 -- particle species: protons, solar activity: min !! * 4 -- particle species: protons, solar activity: max - real(wp) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. + real(wp) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. - real(wp) :: xl !! l value - real(wp) :: bbx + real(wp) :: xl !! l value + real(wp) :: bbx - call me%igrf%igrfc(v,year,xl,bbx) - call me%trm%aep8(e,xl,bbx,imname,flux) + call me%igrf%igrfc(v, year, xl, bbx) + call me%trm%aep8(e, xl, bbx, imname, flux) - end function get_flux_c_ + end function get_flux_c_ !***************************************************************************************** !***************************************************************************************** @@ -181,24 +181,24 @@ end function get_flux_c_ !@note This routine is not efficient at all since it will reload all the ! files every time it is called. - function get_flux_c(v,year,e,imname) result(flux) + function get_flux_c(v, year, e, imname) result(flux) - real(wp),dimension(3),intent(in) :: v - real(wp),intent(in) :: year !! decimal year for which geomagnetic field is to + real(wp), dimension(3), intent(in) :: v + real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) - real(wp),intent(in) :: e !! minimum energy - integer,intent(in) :: imname !! which method to use: + real(wp), intent(in) :: e !! minimum energy + integer, intent(in) :: imname !! which method to use: !! !! * 1 -- particle species: electrons, solar activity: min !! * 2 -- particle species: electrons, solar activity: max !! * 3 -- particle species: protons, solar activity: min !! * 4 -- particle species: protons, solar activity: max - real(wp) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. + real(wp) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1. - type(radbelt_type) :: radbelt + type(radbelt_type) :: radbelt - flux = radbelt%get_flux(v,year,e,imname) + flux = radbelt%get_flux(v, year, e, imname) - end function get_flux_c + end function get_flux_c end module radbelt_module diff --git a/src/shellig.f90 b/src/shellig.f90 index 388d230..232016a 100644 --- a/src/shellig.f90 +++ b/src/shellig.f90 @@ -11,302 +11,302 @@ ! * 5/31/00-DKB- Updated to IGRF-45-00; new coeff.: IGRF00, IGRF00s ! * 3/24/05-DKB- Updated to IGRF-45-10; new coeff.: IGRF05, IGRF05s - module shellig_module +module shellig_module - use radbelt_kinds_module + use radbelt_kinds_module - implicit none + implicit none - private + private - integer,parameter :: filename_len = 14 !! length of the model data file names + integer, parameter :: filename_len = 14 !! length of the model data file names - ! parameters formerly in `gener` common block - real(wp),parameter :: Era = 6371.2_wp !! earth radius for normalization of cartesian coordinates (6371.2 km) - real(wp),parameter :: erequ = 6378.16_wp - real(wp),parameter :: erpol = 6356.775_wp - real(wp),parameter :: Aquad = erequ*erequ !! square of major half axis for + ! parameters formerly in `gener` common block + real(wp), parameter :: Era = 6371.2_wp !! earth radius for normalization of cartesian coordinates (6371.2 km) + real(wp), parameter :: erequ = 6378.16_wp + real(wp), parameter :: erpol = 6356.775_wp + real(wp), parameter :: Aquad = erequ * erequ !! square of major half axis for !! earth ellipsoid as recommended by international !! astronomical union - real(wp),parameter :: Bquad = erpol*erpol !! square of minor half axis for + real(wp), parameter :: Bquad = erpol * erpol !! square of minor half axis for !! earth ellipsoid as recommended by international !! astronomical union - real(wp),parameter :: Umr = atan(1.0_wp)*4.0_wp/180.0_wp !! atan(1.0)*4./180. *umr= + real(wp), parameter :: Umr = atan(1.0_wp) * 4.0_wp / 180.0_wp !! atan(1.0)*4./180. *umr= - real(wp),dimension(3,3),parameter :: u = reshape([ +0.3511737_wp , -0.9148385_wp , -0.1993679_wp , & - +0.9335804_wp , +0.3583680_wp , +0.0000000_wp , & - +0.0714471_wp , -0.1861260_wp , +0.9799247_wp], [3,3]) - integer,parameter :: max_loop_index = 3333 !! used in [[shellg]] for the field line integration loop + real(wp), dimension(3, 3), parameter :: u = reshape([+0.3511737_wp, -0.9148385_wp, -0.1993679_wp, & + +0.9335804_wp, +0.3583680_wp, +0.0000000_wp, & + +0.0714471_wp, -0.1861260_wp, +0.9799247_wp], [3, 3]) + integer, parameter :: max_loop_index = 3333 !! used in [[shellg]] for the field line integration loop - type,public :: shellig_type - private + type, public :: shellig_type + private - character(len=:),allocatable :: igrf_dir !! directory containing the data files + character(len=:), allocatable :: igrf_dir !! directory containing the data files - ! formerly in the `fidb0` common block - real(wp),dimension(3) :: sp = 0.0_wp + ! formerly in the `fidb0` common block + real(wp), dimension(3) :: sp = 0.0_wp - ! formerly in blank common - real(wp),dimension(3) :: xi = 0.0_wp - real(wp),dimension(144) :: h = 0.0_wp !! Field model coefficients adjusted for [[shellg]] + ! formerly in blank common + real(wp), dimension(3) :: xi = 0.0_wp + real(wp), dimension(144) :: h = 0.0_wp !! Field model coefficients adjusted for [[shellg]] - ! formerly in `model` common block - integer :: iyea = 0 !! the int year corresponding to the file `name` that has been read - character(len=:),allocatable :: name !! file name - integer :: nmax = 0 !! maximum order of spherical harmonics - real(wp) :: Time = 0.0_wp !! year (decimal: 1973.5) for which magnetic field is to be calculated - real(wp),dimension(144) :: g = 0.0_wp !! `g(m)` -- normalized field coefficients (see [[feldcof]]) m=nmax*(nmax+2) - integer :: nmax1 = 0 !! saved variables from the file - integer :: nmax2 = 0 !! saved variables from the file - real(wp),dimension(144) :: g_cache = 0.0_wp !! saved `g` from the file + ! formerly in `model` common block + integer :: iyea = 0 !! the int year corresponding to the file `name` that has been read + character(len=:), allocatable :: name !! file name + integer :: nmax = 0 !! maximum order of spherical harmonics + real(wp) :: Time = 0.0_wp !! year (decimal: 1973.5) for which magnetic field is to be calculated + real(wp), dimension(144) :: g = 0.0_wp !! `g(m)` -- normalized field coefficients (see [[feldcof]]) m=nmax*(nmax+2) + integer :: nmax1 = 0 !! saved variables from the file + integer :: nmax2 = 0 !! saved variables from the file + real(wp), dimension(144) :: g_cache = 0.0_wp !! saved `g` from the file - ! formerly saved vars in shellg: - real(wp) :: step = 0.20_wp !! step size for field line tracing - real(wp) :: steq = 0.03_wp !! step size for integration + ! formerly saved vars in shellg: + real(wp) :: step = 0.20_wp !! step size for field line tracing + real(wp) :: steq = 0.03_wp !! step size for integration - ! from feldcof, so we can cache the coefficients - real(wp),dimension(120) :: gh2 = 0.0_wp ! JW : why is this 120 and g is 144 ??? + ! from feldcof, so we can cache the coefficients + real(wp), dimension(120) :: gh2 = 0.0_wp ! JW : why is this 120 and g is 144 ??? - real(wp),dimension(:,:),allocatable :: p !! this was `p(8,100)` in the original code. + real(wp), dimension(:, :), allocatable :: p !! this was `p(8,100)` in the original code. !! used for the field line integration loop. !! changed it to be allocatable since it was !! changed to be p(8,3334). - contains - private + contains + private - procedure,public :: igrf, igrfc + procedure, public :: igrf, igrfc - procedure, public :: feldcof - procedure, public :: feldg, feldc - procedure, public :: shellg, shellc - procedure, public :: findb0 - procedure :: stoer, feldi - procedure,public :: set_data_file_dir, get_data_file_dir - procedure,public :: destroy => destroy_shellig_type + procedure, public :: feldcof + procedure, public :: feldg, feldc + procedure, public :: shellg, shellc + procedure, public :: findb0 + procedure :: stoer, feldi + procedure, public :: set_data_file_dir, get_data_file_dir + procedure, public :: destroy => destroy_shellig_type - end type shellig_type + end type shellig_type - contains +contains !***************************************************************************************** !***************************************************************************************** !> ! Destroy a [[shellig_type]]. - subroutine destroy_shellig_type(me) - class(shellig_type),intent(out) :: me - end subroutine destroy_shellig_type + subroutine destroy_shellig_type(me) + class(shellig_type), intent(out) :: me + end subroutine destroy_shellig_type !***************************************************************************************** !> ! Set the directory containing the data files. - subroutine set_data_file_dir(me,dir) - class(shellig_type),intent(inout) :: me - character(len=*),intent(in) :: dir - me%igrf_dir = trim(dir) - end subroutine set_data_file_dir + subroutine set_data_file_dir(me, dir) + class(shellig_type), intent(inout) :: me + character(len=*), intent(in) :: dir + me%igrf_dir = trim(dir) + end subroutine set_data_file_dir !***************************************************************************************** !> ! Get the directory containing the data files. - function get_data_file_dir(me) result(dir) - class(shellig_type),intent(in) :: me - character(len=:),allocatable :: dir - if (allocated(me%igrf_dir)) then - dir = trim(me%igrf_dir) // '/' - else - dir = 'data/igrf/' ! default - end if - end function get_data_file_dir + function get_data_file_dir(me) result(dir) + class(shellig_type), intent(in) :: me + character(len=:), allocatable :: dir + if (allocated(me%igrf_dir)) then + dir = trim(me%igrf_dir)//'/' + else + dir = 'data/igrf/' ! default + end if + end function get_data_file_dir !***************************************************************************************** !> ! Wrapper for IGRF functions. - subroutine igrf(me,lon,lat,height,year,xl,bbx) + subroutine igrf(me, lon, lat, height, year, xl, bbx) - class(shellig_type),intent(inout) :: me - real(wp),intent(in) :: lon !! geodetic longitude in degrees (east) - real(wp),intent(in) :: lat !! geodetic latitude in degrees (north) - real(wp),intent(in) :: height !! altitude in km above sea level - real(wp),intent(in) :: year !! decimal year for which geomagnetic field is to + class(shellig_type), intent(inout) :: me + real(wp), intent(in) :: lon !! geodetic longitude in degrees (east) + real(wp), intent(in) :: lat !! geodetic latitude in degrees (north) + real(wp), intent(in) :: height !! altitude in km above sea level + real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) - real(wp),intent(out) :: xl !! l-value - real(wp),intent(out) :: bbx !! b_total / b_equatorial ratio - - real(wp) :: bab1 , babs , bdel , bdown , beast , & - beq , bequ , bnorth , dimo , rr0 - integer :: icode - logical :: val - - real(wp),parameter :: stps = 0.05_wp - - ! JW : do we need to reset some or all of these ? - me%sp = 0.0_wp - me%xi = 0.0_wp - me%h = 0.0_wp - me%step = 0.20_wp - me%steq = 0.03_wp - - call me%feldcof(year,dimo) - call me%feldg(lat,lon,height,bnorth,beast,bdown,babs) - call me%shellg(lat,lon,height,dimo,xl,icode,bab1) - - bequ = dimo/(xl*xl*xl) - if ( icode==1 ) then - bdel = 1.0e-3_wp - call me%findb0(stps,bdel,val,beq,rr0) - if ( val ) bequ = beq - endif - bbx = babs/bequ - - end subroutine igrf + real(wp), intent(out) :: xl !! l-value + real(wp), intent(out) :: bbx !! b_total / b_equatorial ratio + + real(wp) :: bab1, babs, bdel, bdown, beast, & + beq, bequ, bnorth, dimo, rr0 + integer :: icode + logical :: val + + real(wp), parameter :: stps = 0.05_wp + + ! JW : do we need to reset some or all of these ? + me%sp = 0.0_wp + me%xi = 0.0_wp + me%h = 0.0_wp + me%step = 0.20_wp + me%steq = 0.03_wp + + call me%feldcof(year, dimo) + call me%feldg(lat, lon, height, bnorth, beast, bdown, babs) + call me%shellg(lat, lon, height, dimo, xl, icode, bab1) + + bequ = dimo / (xl * xl * xl) + if (icode == 1) then + bdel = 1.0e-3_wp + call me%findb0(stps, bdel, val, beq, rr0) + if (val) bequ = beq + end if + bbx = babs / bequ + + end subroutine igrf !***************************************************************************************** !***************************************************************************************** !> ! Alternate version of [[igrf]] for cartesian coordinates. - subroutine igrfc(me,v,year,xl,bbx) + subroutine igrfc(me, v, year, xl, bbx) - class(shellig_type),intent(inout) :: me - real(wp),dimension(3),intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) + class(shellig_type), intent(inout) :: me + real(wp), dimension(3), intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) !! x-axis pointing to equator at 0 longitude !! y-axis pointing to equator at 90 long. !! z-axis pointing to north pole - real(wp),intent(in) :: year !! decimal year for which geomagnetic field is to + real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) - real(wp),intent(out) :: xl !! l-value - real(wp),intent(out) :: bbx !! b_total / b_equatorial ratio - - real(wp) :: bab1 , bdel , beq , bequ , dimo , rr0 - integer :: icode - logical :: val - real(wp),dimension(3) :: b - - real(wp),parameter :: stps = 0.05_wp - - ! JW : do we need to reset some or all of these ? - me%sp = 0.0_wp - me%xi = 0.0_wp - me%h = 0.0_wp - me%step = 0.20_wp - me%steq = 0.03_wp - - call me%feldcof(year,dimo) - call me%feldc(v,b) - call me%shellc(v,dimo,xl,icode,bab1) - - bequ = dimo/(xl*xl*xl) - if ( icode==1 ) then - bdel = 1.0e-3_wp - call me%findb0(stps,bdel,val,beq,rr0) - if ( val ) bequ = beq - endif - bbx = norm2(b)/bequ - - end subroutine igrfc + real(wp), intent(out) :: xl !! l-value + real(wp), intent(out) :: bbx !! b_total / b_equatorial ratio + + real(wp) :: bab1, bdel, beq, bequ, dimo, rr0 + integer :: icode + logical :: val + real(wp), dimension(3) :: b + + real(wp), parameter :: stps = 0.05_wp + + ! JW : do we need to reset some or all of these ? + me%sp = 0.0_wp + me%xi = 0.0_wp + me%h = 0.0_wp + me%step = 0.20_wp + me%steq = 0.03_wp + + call me%feldcof(year, dimo) + call me%feldc(v, b) + call me%shellc(v, dimo, xl, icode, bab1) + + bequ = dimo / (xl * xl * xl) + if (icode == 1) then + bdel = 1.0e-3_wp + call me%findb0(stps, bdel, val, beq, rr0) + if (val) bequ = beq + end if + bbx = norm2(b) / bequ + + end subroutine igrfc !***************************************************************************************** !***************************************************************************************** !> - subroutine findb0(me,stps,bdel,value,bequ,rr0) - - class(shellig_type),intent(inout) :: me - real(wp),intent(in) :: stps - real(wp),intent(inout) :: bdel - real(wp),intent(out) :: bequ - logical,intent(out) :: value - real(wp),intent(out) :: rr0 - - real(wp) :: b , bdelta , bmin , bold , bq1 , & - bq2 , bq3 , p(8,4) , r1 , r2 , r3 , & - rold , step , step12 , zz - integer :: i , irun , j , n - - step=stps - irun=0 - rold = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings - - main : do - irun=irun+1 - if (irun>5) then - value=.false. - exit main - endif - ! first three points - p(1,2)=me%sp(1) - p(2,2)=me%sp(2) - p(3,2)=me%sp(3) - step=-sign(step,p(3,2)) - call me%stoer(p(1,2),bq2,r2) - p(1,3)=p(1,2)+0.5_wp*step*p(4,2) - p(2,3)=p(2,2)+0.5_wp*step*p(5,2) - p(3,3)=p(3,2)+0.5_wp*step - call me%stoer(p(1,3),bq3,r3) - p(1,1)=p(1,2)-step*(2.0_wp*p(4,2)-p(4,3)) - p(2,1)=p(2,2)-step*(2.0_wp*p(5,2)-p(5,3)) - p(3,1)=p(3,2)-step - call me%stoer(p(1,1),bq1,r1) - p(1,3)=p(1,2)+step*(20.0_wp*p(4,3)-3.*p(4,2)+p(4,1))/18.0_wp - p(2,3)=p(2,2)+step*(20.0_wp*p(5,3)-3.*p(5,2)+p(5,1))/18.0_wp - p(3,3)=p(3,2)+step - call me%stoer(p(1,3),bq3,r3) - ! invert sense if required - if (bq3>bq1) then - step=-step - r3=r1 - bq3=bq1 - do i=1,5 - zz=p(i,1) - p(i,1)=p(i,3) - p(i,3)=zz - end do - end if - ! initialization - step12=step/12.0_wp - value=.true. - bmin=1.0e4_wp - bold=1.0e4_wp - ! corrector (field line tracing) - n=0 - corrector : do - p(1,3)=p(1,2)+step12*(5.0_wp*p(4,3)+8.0_wp*p(4,2)-p(4,1)) - n=n+1 - p(2,3)=p(2,2)+step12*(5.0_wp*p(5,3)+8.0_wp*p(5,2)-p(5,1)) - ! predictor (field line tracing) - p(1,4)=p(1,3)+step12*(23.0_wp*p(4,3)-16.0_wp*p(4,2)+5.0_wp*p(4,1)) - p(2,4)=p(2,3)+step12*(23.0_wp*p(5,3)-16.0_wp*p(5,2)+5.0_wp*p(5,1)) - p(3,4)=p(3,3)+step - call me%stoer(p(1,4),bq3,r3) - do j=1,3 - do i=1,8 - p(i,j)=p(i,j+1) + subroutine findb0(me, stps, bdel, value, bequ, rr0) + + class(shellig_type), intent(inout) :: me + real(wp), intent(in) :: stps + real(wp), intent(inout) :: bdel + real(wp), intent(out) :: bequ + logical, intent(out) :: value + real(wp), intent(out) :: rr0 + + real(wp) :: b, bdelta, bmin, bold, bq1, & + bq2, bq3, p(8, 4), r1, r2, r3, & + rold, step, step12, zz + integer :: i, irun, j, n + + step = stps + irun = 0 + rold = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings + + main: do + irun = irun + 1 + if (irun > 5) then + value = .false. + exit main + end if + ! first three points + p(1, 2) = me%sp(1) + p(2, 2) = me%sp(2) + p(3, 2) = me%sp(3) + step = -sign(step, p(3, 2)) + call me%stoer(p(1, 2), bq2, r2) + p(1, 3) = p(1, 2) + 0.5_wp * step * p(4, 2) + p(2, 3) = p(2, 2) + 0.5_wp * step * p(5, 2) + p(3, 3) = p(3, 2) + 0.5_wp * step + call me%stoer(p(1, 3), bq3, r3) + p(1, 1) = p(1, 2) - step * (2.0_wp * p(4, 2) - p(4, 3)) + p(2, 1) = p(2, 2) - step * (2.0_wp * p(5, 2) - p(5, 3)) + p(3, 1) = p(3, 2) - step + call me%stoer(p(1, 1), bq1, r1) + p(1, 3) = p(1, 2) + step * (20.0_wp * p(4, 3) - 3.*p(4, 2) + p(4, 1)) / 18.0_wp + p(2, 3) = p(2, 2) + step * (20.0_wp * p(5, 3) - 3.*p(5, 2) + p(5, 1)) / 18.0_wp + p(3, 3) = p(3, 2) + step + call me%stoer(p(1, 3), bq3, r3) + ! invert sense if required + if (bq3 > bq1) then + step = -step + r3 = r1 + bq3 = bq1 + do i = 1, 5 + zz = p(i, 1) + p(i, 1) = p(i, 3) + p(i, 3) = zz end do - end do - b=sqrt(bq3) - if (bbold) exit corrector - bold=b - rold=1.0_wp/r3 - me%sp(1)=p(1,4) - me%sp(2)=p(2,4) - me%sp(3)=p(3,4) - end do corrector - if (bold/=bmin) value=.false. - bdelta=(b-bold)/bold - if (bdelta<=bdel) exit main - step=step/10.0_wp - end do main - - rr0=rold - bequ=bold - bdel=bdelta - - end subroutine findb0 + end if + ! initialization + step12 = step / 12.0_wp + value = .true. + bmin = 1.0e4_wp + bold = 1.0e4_wp + ! corrector (field line tracing) + n = 0 + corrector: do + p(1, 3) = p(1, 2) + step12 * (5.0_wp * p(4, 3) + 8.0_wp * p(4, 2) - p(4, 1)) + n = n + 1 + p(2, 3) = p(2, 2) + step12 * (5.0_wp * p(5, 3) + 8.0_wp * p(5, 2) - p(5, 1)) + ! predictor (field line tracing) + p(1, 4) = p(1, 3) + step12 * (23.0_wp * p(4, 3) - 16.0_wp * p(4, 2) + 5.0_wp * p(4, 1)) + p(2, 4) = p(2, 3) + step12 * (23.0_wp * p(5, 3) - 16.0_wp * p(5, 2) + 5.0_wp * p(5, 1)) + p(3, 4) = p(3, 3) + step + call me%stoer(p(1, 4), bq3, r3) + do j = 1, 3 + do i = 1, 8 + p(i, j) = p(i, j + 1) + end do + end do + b = sqrt(bq3) + if (b < bmin) bmin = b + if (b > bold) exit corrector + bold = b + rold = 1.0_wp / r3 + me%sp(1) = p(1, 4) + me%sp(2) = p(2, 4) + me%sp(3) = p(3, 4) + end do corrector + if (bold /= bmin) value = .false. + bdelta = (b - bold) / bold + if (bdelta <= bdel) exit main + step = step / 10.0_wp + end do main + + rr0 = rold + bequ = bold + bdel = bdelta + + end subroutine findb0 !***************************************************************************************** !> @@ -315,26 +315,26 @@ end subroutine findb0 !@note In the original code, this was an ENTRY point in [[shellg]] and didn't ! include all the outputs. - subroutine shellc(me,v,dimo,fl,icode,b0) + subroutine shellc(me, v, dimo, fl, icode, b0) - class(shellig_type),intent(inout) :: me - real(wp),dimension(3),intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) + class(shellig_type), intent(inout) :: me + real(wp), dimension(3), intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) !! * x-axis pointing to equator at 0 longitude !! * y-axis pointing to equator at 90 long. !! * z-axis pointing to north pole - real(wp),intent(in) :: dimo !! dipol moment in gauss (normalized to earth radius) - real(wp),intent(out) :: fl !! l-value - integer,intent(out) :: icode !! * =1 normal completion + real(wp), intent(in) :: dimo !! dipol moment in gauss (normalized to earth radius) + real(wp), intent(out) :: fl !! l-value + integer, intent(out) :: icode !! * =1 normal completion !! * =2 unphysical conjugate point (fl meaningless) !! * =3 shell parameter greater than limit up to !! which accurate calculation is required; !! approximation is used. - real(wp),intent(out) :: b0 !! magnetic field strength in gauss - real(wp) :: glat,glon,alt !! not used + real(wp), intent(out) :: b0 !! magnetic field strength in gauss + real(wp) :: glat, glon, alt !! not used - call me%shellg(glat,glon,alt,dimo,fl,icode,b0,v) + call me%shellg(glat, glon, alt, dimo, fl, icode, b0, v) - end subroutine shellc + end subroutine shellc !***************************************************************************************** !> @@ -351,21 +351,21 @@ end subroutine shellc ! - USING CORRECT DIPOL MOMENT I.E.,DIFFERENT COMMON/MODEL/ ! - USING IGRF EARTH MAGNETIC FIELD MODELS FROM 1945 TO 1990 - subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v) + subroutine shellg(me, glat, glon, alt, dimo, fl, icode, b0, v) - class(shellig_type),intent(inout) :: me - real(wp),intent(in) :: glat !! geodetic latitude in degrees (north) - real(wp),intent(in) :: glon !! geodetic longitude in degrees (east) - real(wp),intent(in) :: alt !! altitude in km above sea level - real(wp),intent(in) :: dimo !! dipol moment in gauss (normalized to earth radius) - real(wp),intent(out) :: fl !! l-value - integer,intent(out) :: icode !! * =1 normal completion + class(shellig_type), intent(inout) :: me + real(wp), intent(in) :: glat !! geodetic latitude in degrees (north) + real(wp), intent(in) :: glon !! geodetic longitude in degrees (east) + real(wp), intent(in) :: alt !! altitude in km above sea level + real(wp), intent(in) :: dimo !! dipol moment in gauss (normalized to earth radius) + real(wp), intent(out) :: fl !! l-value + integer, intent(out) :: icode !! * =1 normal completion !! * =2 unphysical conjugate point (fl meaningless) !! * =3 shell parameter greater than limit up to !! which accurate calculation is required; !! approximation is used. - real(wp),intent(out) :: b0 !! magnetic field strength in gauss - real(wp),dimension(3),intent(in),optional :: v !! cartesian coordinates in earth radii (6371.2 km) + real(wp), intent(out) :: b0 !! magnetic field strength in gauss + real(wp), dimension(3), intent(in), optional :: v !! cartesian coordinates in earth radii (6371.2 km) !! !! * x-axis pointing to equator at 0 longitude !! * y-axis pointing to equator at 90 long. @@ -374,251 +374,251 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v) !! If this argument is present, it is used !! instead of glat,glon,alt. See [[shellc]]. - real(wp) :: arg1 , arg2 , bequ , bq1 , bq2 , bq3 , c0 , c1 , c2 , c3 , & - d0 , d1 , d2, dimob0 , e0 , e1 , e2 , ff , fi , gg , & - hli , oradik , oterm , r , r1 , r2 , r3 , r3h , radik , & - rq , step12 , step2 , stp , t , term , xx , z , zq , zz - integer :: i , iequ , n - - real(wp),parameter :: rmin = 0.05_wp !! boundaries for identification of `icode=2 and 3` - real(wp),parameter :: rmax = 1.01_wp !! boundaries for identification of `icode=2 and 3` - - if (.not. allocated(me%p)) allocate(me%p(8,max_loop_index+1)) ! because `p(:,n+1)` in the loop - - bequ = 1.0e10_wp - - if (present(v)) then - me%xi(1) = v(1) - me%xi(2) = v(2) - me%xi(3) = v(3) - else - me%xi = geo_to_cart(glat,glon,alt) - end if - - associate (p => me%p) - - ! convert to dipol-oriented co-ordinates - rq = 1.0_wp/(me%xi(1)*me%xi(1)+me%xi(2)*me%xi(2)+me%xi(3)*me%xi(3)) - r3h = sqrt(rq*sqrt(rq)) - p(1,2) = (me%xi(1)*u(1,1)+me%xi(2)*u(2,1)+me%xi(3)*u(3,1))*r3h - p(2,2) = (me%xi(1)*u(1,2)+me%xi(2)*u(2,2))*r3h - p(3,2) = (me%xi(1)*u(1,3)+me%xi(2)*u(2,3)+me%xi(3)*u(3,3))*rq - ! first three points of field line - me%step = -sign(me%step,p(3,2)) - call me%stoer(p(1,2),bq2,r2) - b0 = sqrt(bq2) - p(1,3) = p(1,2) + 0.5_wp*me%step*p(4,2) - p(2,3) = p(2,2) + 0.5_wp*me%step*p(5,2) - p(3,3) = p(3,2) + 0.5_wp*me%step - call me%stoer(p(1,3),bq3,r3) - p(1,1) = p(1,2) - me%step*(2.0_wp*p(4,2)-p(4,3)) - p(2,1) = p(2,2) - me%step*(2.0_wp*p(5,2)-p(5,3)) - p(3,1) = p(3,2) - me%step - call me%stoer(p(1,1),bq1,r1) - p(1,3) = p(1,2) + me%step*(20.0_wp*p(4,3)-3.*p(4,2)+p(4,1))/18.0_wp - p(2,3) = p(2,2) + me%step*(20.0_wp*p(5,3)-3.*p(5,2)+p(5,1))/18.0_wp - p(3,3) = p(3,2) + me%step - call me%stoer(p(1,3),bq3,r3) - ! invert sense if required - if ( bq3>bq1 ) then - me%step = -me%step - r3 = r1 - bq3 = bq1 - do i = 1 , 7 - zz = p(i,1) - p(i,1) = p(i,3) - p(i,3) = zz - enddo - endif - ! search for lowest magnetic field strength - if ( bq11.0_wp ) then - ! predictor (field line tracing) - p(1,n+1) = p(1,n) + step12*(23.0_wp*p(4,n)-16.0_wp*p(4,n-1)+5.0_wp*p(4,n-2)) - p(2,n+1) = p(2,n) + step12*(23.0_wp*p(5,n)-16.0_wp*p(5,n-1)+5.0_wp*p(5,n-2)) - p(3,n+1) = p(3,n) + me%step - call me%stoer(p(1,n+1),bq3,r3) - ! search for lowest magnetic field strength - if ( bq3 me%p) + + ! convert to dipol-oriented co-ordinates + rq = 1.0_wp / (me%xi(1) * me%xi(1) + me%xi(2) * me%xi(2) + me%xi(3) * me%xi(3)) + r3h = sqrt(rq * sqrt(rq)) + p(1, 2) = (me%xi(1) * u(1, 1) + me%xi(2) * u(2, 1) + me%xi(3) * u(3, 1)) * r3h + p(2, 2) = (me%xi(1) * u(1, 2) + me%xi(2) * u(2, 2)) * r3h + p(3, 2) = (me%xi(1) * u(1, 3) + me%xi(2) * u(2, 3) + me%xi(3) * u(3, 3)) * rq + ! first three points of field line + me%step = -sign(me%step, p(3, 2)) + call me%stoer(p(1, 2), bq2, r2) + b0 = sqrt(bq2) + p(1, 3) = p(1, 2) + 0.5_wp * me%step * p(4, 2) + p(2, 3) = p(2, 2) + 0.5_wp * me%step * p(5, 2) + p(3, 3) = p(3, 2) + 0.5_wp * me%step + call me%stoer(p(1, 3), bq3, r3) + p(1, 1) = p(1, 2) - me%step * (2.0_wp * p(4, 2) - p(4, 3)) + p(2, 1) = p(2, 2) - me%step * (2.0_wp * p(5, 2) - p(5, 3)) + p(3, 1) = p(3, 2) - me%step + call me%stoer(p(1, 1), bq1, r1) + p(1, 3) = p(1, 2) + me%step * (20.0_wp * p(4, 3) - 3.*p(4, 2) + p(4, 1)) / 18.0_wp + p(2, 3) = p(2, 2) + me%step * (20.0_wp * p(5, 3) - 3.*p(5, 2) + p(5, 1)) / 18.0_wp + p(3, 3) = p(3, 2) + me%step + call me%stoer(p(1, 3), bq3, r3) + ! invert sense if required + if (bq3 > bq1) then + me%step = -me%step + r3 = r1 + bq3 = bq1 + do i = 1, 7 + zz = p(i, 1) + p(i, 1) = p(i, 3) + p(i, 3) = zz + end do + end if + ! search for lowest magnetic field strength + if (bq1 < bequ) then + bequ = bq1 + iequ = 1 + end if + if (bq2 < bequ) then + bequ = bq2 + iequ = 2 + end if + if (bq3 < bequ) then + bequ = bq3 + iequ = 3 + end if + ! initialization of integration loops + step12 = me%step / 12.0_wp + step2 = me%step + me%step + me%steq = sign(me%steq, me%step) + fi = 0.0_wp + icode = 1 + oradik = 0.0_wp + oterm = 0.0_wp + stp = r2 * me%steq + z = p(3, 2) + stp + stp = stp / 0.75_wp + p(8, 1) = step2 * (p(1, 1) * p(4, 1) + p(2, 1) * p(5, 1)) + p(8, 2) = step2 * (p(1, 2) * p(4, 2) + p(2, 2) * p(5, 2)) + ! main loop (field line tracing) + main: do n = 3, max_loop_index + ! corrector (field line tracing) + p(1, n) = p(1, n - 1) + step12 * (5.0_wp * p(4, n) + 8.0_wp * p(4, n - 1) - p(4, n - 2)) + p(2, n) = p(2, n - 1) + step12 * (5.0_wp * p(5, n) + 8.0_wp * p(5, n - 1) - p(5, n - 2)) + ! prepare expansion coefficients for interpolation + ! of slowly varying quantities + p(8, n) = step2 * (p(1, n) * p(4, n) + p(2, n) * p(5, n)) + c0 = p(1, n - 1)**2 + p(2, n - 1)**2 + c1 = p(8, n - 1) + c2 = (p(8, n) - p(8, n - 2)) * 0.25_wp + c3 = (p(8, n) + p(8, n - 2) - c1 - c1) / 6.0_wp + d0 = p(6, n - 1) + d1 = (p(6, n) - p(6, n - 2)) * 0.5_wp + d2 = (p(6, n) + p(6, n - 2) - d0 - d0) * 0.5_wp + e0 = p(7, n - 1) + e1 = (p(7, n) - p(7, n - 2)) * 0.5_wp + e2 = (p(7, n) + p(7, n - 2) - e0 - e0) * 0.5_wp + inner: do + ! inner loop (for quadrature) + t = (z - p(3, n - 1)) / me%step + if (t > 1.0_wp) then + ! predictor (field line tracing) + p(1, n + 1) = p(1, n) + step12 * (23.0_wp * p(4, n) - 16.0_wp * p(4, n - 1) + 5.0_wp * p(4, n - 2)) + p(2, n + 1) = p(2, n) + step12 * (23.0_wp * p(5, n) - 16.0_wp * p(5, n - 1) + 5.0_wp * p(5, n - 2)) + p(3, n + 1) = p(3, n) + me%step + call me%stoer(p(1, n + 1), bq3, r3) + ! search for lowest magnetic field strength + if (bq3 < bequ) then + iequ = n + 1 + bequ = bq3 + end if + exit inner + else + hli = 0.5_wp * (((c3 * t + c2) * t + c1) * t + c0) + zq = z * z + r = hli + sqrt(hli * hli + zq) + if (r <= rmin) then + ! approximation for high values of l. + icode = 3 + t = -p(3, n - 1) / me%step + fl = 1.0_wp / (abs(((c3 * t + c2) * t + c1) * t + c0) + 1.0e-15_wp) + return + end if + rq = r * r + ff = sqrt(1.0_wp + 3.0_wp * zq / rq) + radik = b0 - ((d2 * t + d1) * t + d0) * r * rq * ff + if (r > rmax) then + icode = 2 + radik = radik - 12.0_wp * (r - rmax)**2 + end if + if (radik + radik <= oradik) exit main + term = sqrt(radik) * ff * ((e2 * t + e1) * t + e0) / (rq + zq) + fi = fi + stp * (oterm + term) + oradik = radik + oterm = term + stp = r * me%steq + z = z + stp + end if + end do inner + end do main + if (iequ < 2) iequ = 2 + me%sp(1) = p(1, iequ - 1) + me%sp(2) = p(2, iequ - 1) + me%sp(3) = p(3, iequ - 1) + if (oradik >= 1.0e-15_wp) fi = fi + stp / 0.75_wp * oterm * oradik / (oradik - radik) + + ! the minimal allowable value of fi was changed from 1e-15 to 1e-12, + ! because 1e-38 is the minimal allowable arg. for alog in our envir. + ! d. bilitza, nov 87. + fi = 0.5_wp * abs(fi) / sqrt(b0) + 1.0e-12_wp + + ! compute l from b and i. same as carmel in invar. + ! correct dipole moment is used here. d. bilitza, nov 87. + dimob0 = dimo / b0 + arg1 = log(fi) + arg2 = log(dimob0) + ! arg = fi*fi*fi/dimob0 + ! if(abs(arg)>88.0_wp) arg=88.0_wp + xx = 3 * arg1 - arg2 + if (xx > 23.0_wp) then + gg = xx - 3.0460681_wp + elseif (xx > 11.7_wp) then + gg = (((((2.8212095e-8_wp * xx - 3.8049276e-6_wp) * xx + & + 2.170224e-4_wp) * xx - 6.7310339e-3_wp) * xx + & + 1.2038224e-1_wp) * xx - 1.8461796e-1_wp) * xx + 2.0007187_wp + elseif (xx > +3.0_wp) then + gg = ((((((((6.3271665e-10_wp * xx - 3.958306e-8_wp) * xx + & + 9.9766148e-07_wp) * xx - 1.2531932e-5_wp) * xx + & + 7.9451313e-5_wp) * xx - 3.2077032e-4_wp) * xx + & + 2.1680398e-3_wp) * xx + 1.2817956e-2_wp) * xx + & + 4.3510529e-1_wp) * xx + 6.222355e-1_wp + elseif (xx > -3.0_wp) then + gg = ((((((((2.6047023e-10_wp * xx + 2.3028767e-9_wp) * xx - & + 2.1997983e-8_wp) * xx - 5.3977642e-7_wp) * xx - & + 3.3408822e-6_wp) * xx + 3.8379917e-5_wp) * xx + & + 1.1784234e-3_wp) * xx + 1.4492441e-2_wp) * xx + & + 4.3352788e-1_wp) * xx + 6.228644e-1_wp + elseif (xx > -22.0_wp) then + gg = ((((((((-8.1537735e-14_wp * xx + 8.3232531e-13_wp) * xx + & + 1.0066362e-9_wp) * xx + 8.1048663e-8_wp) * xx + & + 3.2916354e-6_wp) * xx + 8.2711096e-5_wp) * xx + & + 1.3714667e-3_wp) * xx + 1.5017245e-2_wp) * xx + & + 4.3432642e-1_wp) * xx + 6.2337691e-1_wp else - hli = 0.5_wp*(((c3*t+c2)*t+c1)*t+c0) - zq = z*z - r = hli + sqrt(hli*hli+zq) - if ( r<=rmin ) then - ! approximation for high values of l. - icode = 3 - t = -p(3,n-1)/me%step - fl = 1.0_wp/(abs(((c3*t+c2)*t+c1)*t+c0)+1.0e-15_wp) - return - endif - rq = r*r - ff = sqrt(1.0_wp+3.0_wp*zq/rq) - radik = b0 - ((d2*t+d1)*t+d0)*r*rq*ff - if ( r>rmax ) then - icode = 2 - radik = radik - 12.0_wp*(r-rmax)**2 - endif - if ( radik+radik<=oradik ) exit main - term = sqrt(radik)*ff*((e2*t+e1)*t+e0)/(rq+zq) - fi = fi + stp*(oterm+term) - oradik = radik - oterm = term - stp = r*me%steq - z = z + stp - endif - enddo inner - enddo main - if ( iequ<2 ) iequ = 2 - me%sp(1) = p(1,iequ-1) - me%sp(2) = p(2,iequ-1) - me%sp(3) = p(3,iequ-1) - if ( oradik>=1.0e-15_wp ) fi = fi + stp/0.75_wp*oterm*oradik/(oradik-radik) - - ! the minimal allowable value of fi was changed from 1e-15 to 1e-12, - ! because 1e-38 is the minimal allowable arg. for alog in our envir. - ! d. bilitza, nov 87. - fi = 0.5_wp*abs(fi)/sqrt(b0) + 1.0e-12_wp - - ! compute l from b and i. same as carmel in invar. - ! correct dipole moment is used here. d. bilitza, nov 87. - dimob0 = dimo/b0 - arg1 = log(fi) - arg2 = log(dimob0) - ! arg = fi*fi*fi/dimob0 - ! if(abs(arg)>88.0_wp) arg=88.0_wp - xx = 3*arg1 - arg2 - if ( xx>23.0_wp ) then - gg = xx - 3.0460681_wp - elseif ( xx>11.7_wp ) then - gg = (((((2.8212095e-8_wp*xx-3.8049276e-6_wp)*xx+& - 2.170224e-4_wp)*xx-6.7310339e-3_wp)*xx+& - 1.2038224e-1_wp)*xx-1.8461796e-1_wp)*xx + 2.0007187_wp - elseif ( xx>+3.0_wp ) then - gg = ((((((((6.3271665e-10_wp*xx-3.958306e-8_wp)*xx+& - 9.9766148e-07_wp)*xx-1.2531932e-5_wp)*xx+& - 7.9451313e-5_wp)*xx-3.2077032e-4_wp)*xx+& - 2.1680398e-3_wp)*xx+1.2817956e-2_wp)*xx+& - 4.3510529e-1_wp)*xx + 6.222355e-1_wp - elseif ( xx>-3.0_wp ) then - gg = ((((((((2.6047023e-10_wp*xx+2.3028767e-9_wp)*xx-& - 2.1997983e-8_wp)*xx-5.3977642e-7_wp)*xx-& - 3.3408822e-6_wp)*xx+3.8379917e-5_wp)*xx+& - 1.1784234e-3_wp)*xx+1.4492441e-2_wp)*xx+& - 4.3352788e-1_wp)*xx + 6.228644e-1_wp - elseif ( xx>-22.0_wp ) then - gg = ((((((((-8.1537735e-14_wp*xx+8.3232531e-13_wp)*xx+& - 1.0066362e-9_wp)*xx+8.1048663e-8_wp)*xx+& - 3.2916354e-6_wp)*xx+8.2711096e-5_wp)*xx+& - 1.3714667e-3_wp)*xx+1.5017245e-2_wp)*xx+& - 4.3432642e-1_wp)*xx + 6.2337691e-1_wp - else - gg = 3.33338e-1_wp*xx + 3.0062102e-1_wp - endif - fl = exp(log((1.0_wp+exp(gg))*dimob0)/3.0_wp) - - end associate - -end subroutine shellg + gg = 3.33338e-1_wp * xx + 3.0062102e-1_wp + end if + fl = exp(log((1.0_wp + exp(gg)) * dimob0) / 3.0_wp) + + end associate + + end subroutine shellg !***************************************************************************************** !> ! subroutine used for field line tracing in [[shellg]]. ! calls entry point [[feldi]] in geomagnetic field subroutine [[feldg]] -subroutine stoer(me,p,bq,r) - - class(shellig_type),intent(inout) :: me - real(wp),dimension(7),intent(inout) :: p - real(wp),intent(out) :: bq - real(wp),intent(out) :: r - - real(wp) :: dr , dsq , dx , dxm , dy , dym , dz , & - dzm , fli , q , rq , wr , xm , ym , zm - - ! xm,ym,zm are geomagnetic cartesian inverse co-ordinates - zm = P(3) - fli = P(1)*P(1) + P(2)*P(2) + 1.0e-15_wp - R = 0.5_wp*(fli+sqrt(fli*fli+(zm+zm)**2)) - rq = R*R - wr = sqrt(R) - xm = P(1)*wr - ym = P(2)*wr - ! transform to geographic co-ordinate system - me%Xi(1) = xm*u(1,1) + ym*u(1,2) + zm*u(1,3) - me%Xi(2) = xm*u(2,1) + ym*u(2,2) + zm*u(2,3) - me%Xi(3) = xm*u(3,1) + zm*u(3,3) - ! compute derivatives - ! Changed from CALL FELDI(XI,H); XI, H are in COMMON block; results - ! are the same; dkb Feb 1998. - ! JW : feb 2024 : xi, h now class variables. - call me%feldi() - q = me%H(1)/rq - dx = me%H(3) + me%H(3) + q*me%Xi(1) - dy = me%H(4) + me%H(4) + q*me%Xi(2) - dz = me%H(2) + me%H(2) + q*me%Xi(3) - ! transform back to geomagnetic co-ordinate system - dxm = u(1,1)*dx + u(2,1)*dy + u(3,1)*dz - dym = u(1,2)*dx + u(2,2)*dy - dzm = u(1,3)*dx + u(2,3)*dy + u(3,3)*dz - dr = (xm*dxm+ym*dym+zm*dzm)/R - ! form slowly varying expressions - P(4) = (wr*dxm-0.5_wp*P(1)*dr)/(R*dzm) - P(5) = (wr*dym-0.5_wp*P(2)*dr)/(R*dzm) - dsq = rq*(dxm*dxm+dym*dym+dzm*dzm) - Bq = dsq*rq*rq - P(6) = sqrt(dsq/(rq+3.0_wp*zm*zm)) - P(7) = P(6)*(rq+zm*zm)/(rq*dzm) - -end subroutine stoer + subroutine stoer(me, p, bq, r) + + class(shellig_type), intent(inout) :: me + real(wp), dimension(7), intent(inout) :: p + real(wp), intent(out) :: bq + real(wp), intent(out) :: r + + real(wp) :: dr, dsq, dx, dxm, dy, dym, dz, & + dzm, fli, q, rq, wr, xm, ym, zm + + ! xm,ym,zm are geomagnetic cartesian inverse co-ordinates + zm = P(3) + fli = P(1) * P(1) + P(2) * P(2) + 1.0e-15_wp + R = 0.5_wp * (fli + sqrt(fli * fli + (zm + zm)**2)) + rq = R * R + wr = sqrt(R) + xm = P(1) * wr + ym = P(2) * wr + ! transform to geographic co-ordinate system + me%Xi(1) = xm * u(1, 1) + ym * u(1, 2) + zm * u(1, 3) + me%Xi(2) = xm * u(2, 1) + ym * u(2, 2) + zm * u(2, 3) + me%Xi(3) = xm * u(3, 1) + zm * u(3, 3) + ! compute derivatives + ! Changed from CALL FELDI(XI,H); XI, H are in COMMON block; results + ! are the same; dkb Feb 1998. + ! JW : feb 2024 : xi, h now class variables. + call me%feldi() + q = me%H(1) / rq + dx = me%H(3) + me%H(3) + q * me%Xi(1) + dy = me%H(4) + me%H(4) + q * me%Xi(2) + dz = me%H(2) + me%H(2) + q * me%Xi(3) + ! transform back to geomagnetic co-ordinate system + dxm = u(1, 1) * dx + u(2, 1) * dy + u(3, 1) * dz + dym = u(1, 2) * dx + u(2, 2) * dy + dzm = u(1, 3) * dx + u(2, 3) * dy + u(3, 3) * dz + dr = (xm * dxm + ym * dym + zm * dzm) / R + ! form slowly varying expressions + P(4) = (wr * dxm - 0.5_wp * P(1) * dr) / (R * dzm) + P(5) = (wr * dym - 0.5_wp * P(2) * dr) / (R * dzm) + dsq = rq * (dxm * dxm + dym * dym + dzm * dzm) + Bq = dsq * rq * rq + P(6) = sqrt(dsq / (rq + 3.0_wp * zm * zm)) + P(7) = P(6) * (rq + zm * zm) / (rq * dzm) + + end subroutine stoer !***************************************************************************************** !> @@ -636,200 +636,200 @@ end subroutine stoer !@note In the original code, [[feldc] and [[feldi]] were ! ENTRY points to this routine - subroutine feldg(me,glat,glon,alt,bnorth,beast,bdown,babs) + subroutine feldg(me, glat, glon, alt, bnorth, beast, bdown, babs) - class(shellig_type),intent(inout) :: me - real(wp),intent(in) :: glat !! geodetic latitude in degrees (north) - real(wp),intent(in) :: glon !! geodetic longitude in degrees (east) - real(wp),intent(in) :: alt !! altitude in km above sea level - real(wp),intent(out) :: bnorth, beast, bdown !! components of the field with respect + class(shellig_type), intent(inout) :: me + real(wp), intent(in) :: glat !! geodetic latitude in degrees (north) + real(wp), intent(in) :: glon !! geodetic longitude in degrees (east) + real(wp), intent(in) :: alt !! altitude in km above sea level + real(wp), intent(out) :: bnorth, beast, bdown !! components of the field with respect !! to the local geodetic coordinate system, with axis !! pointing in the tangential plane to the north, east !! and downward. - real(wp),intent(out) :: Babs !! magnetic field strength in gauss - - real(wp) :: brho , bxxx , byyy , bzzz , cp , ct , d , f , rho , & - rlat , rlon , rq , s , sp , st , t , & - x , xxx , y , yyy , z , zzz - integer :: i , ih , ihmax , il , imax , k , last , m - - ! same calculation as geo_to_cart, but not used here - ! because the intermediate variables are also used below. - rlat = glat*umr - ct = sin(rlat) - st = cos(rlat) - d = sqrt(aquad-(aquad-bquad)*ct*ct) - rlon = glon*umr - cp = cos(rlon) - sp = sin(rlon) - zzz = (alt+bquad/d)*ct/era - rho = (alt+aquad/d)*st/era - xxx = rho*cp - yyy = rho*sp - - rq = 1.0_wp/(xxx*xxx+yyy*yyy+zzz*zzz) - me%xi = [xxx,yyy,zzz] * rq - - ihmax=me%nmax*me%nmax+1 - last=ihmax+me%nmax+me%nmax - imax=me%nmax+me%nmax-1 - do i=ihmax,last - me%h(i)=me%g(i) - end do - do k=1,3,2 - i=imax - ih=ihmax - do - il=ih-i - f=2.0_wp/real(i-k+2, wp) - x=me%xi(1)*f - y=me%xi(2)*f - z=me%xi(3)*(f+f) - i=i-2 - if ((i-1)>=0) then - if ((i-1)>0) then - do m=3,i,2 - me%h(il+m+1)=me%g(il+m+1)+z*me%h(ih+m+1)+x*(me%h(ih+m+3)-& - me%h(ih+m-1))-y*(me%h(ih+m+2)+me%h(ih+m-2)) - me%h(il+m)=me%g(il+m)+z*me%h(ih+m)+x*(me%h(ih+m+2)-& - me%h(ih+m-2))+y*(me%h(ih+m+3)+me%h(ih+m-1)) - end do - end if - me%h(il+2)=me%g(il+2)+z*me%h(ih+2)+x*me%h(ih+4)-y*(me%h(ih+3)+me%h(ih)) - me%h(il+1)=me%g(il+1)+z*me%h(ih+1)+y*me%h(ih+4)+x*(me%h(ih+3)-me%h(ih)) - end if - me%h(il)=me%g(il)+z*me%h(ih)+2.0_wp*(x*me%h(ih+1)+y*me%h(ih+2)) - ih=il - if (i= 0) then + if ((i - 1) > 0) then + do m = 3, i, 2 + me%h(il + m + 1) = me%g(il + m + 1) + z * me%h(ih + m + 1) + x * (me%h(ih + m + 3) - & + me%h(ih + m - 1)) - y * (me%h(ih + m + 2) + me%h(ih + m - 2)) + me%h(il + m) = me%g(il + m) + z * me%h(ih + m) + x * (me%h(ih + m + 2) - & + me%h(ih + m - 2)) + y * (me%h(ih + m + 3) + me%h(ih + m - 1)) + end do + end if + me%h(il + 2) = me%g(il + 2) + z * me%h(ih + 2) + x * me%h(ih + 4) - y * (me%h(ih + 3) + me%h(ih)) + me%h(il + 1) = me%g(il + 1) + z * me%h(ih + 1) + y * me%h(ih + 4) + x * (me%h(ih + 3) - me%h(ih)) + end if + me%h(il) = me%g(il) + z * me%h(ih) + 2.0_wp * (x * me%h(ih + 1) + y * me%h(ih + 2)) + ih = il + if (i < k) exit + end do + end do + + s = 0.5_wp * me%h(1) + 2.0_wp * (me%h(2) * me%xi(3) + me%h(3) * me%xi(1) + me%h(4) * me%xi(2)) + t = (rq + rq) * sqrt(rq) + bxxx = t * (me%h(3) - s * xxx) + byyy = t * (me%h(4) - s * yyy) + bzzz = t * (me%h(2) - s * zzz) + + babs = sqrt(bxxx * bxxx + byyy * byyy + bzzz * bzzz) + beast = byyy * cp - bxxx * sp + brho = byyy * sp + bxxx * cp + bnorth = bzzz * st - brho * ct + bdown = -bzzz * ct - brho * st + + end subroutine feldg !***************************************************************************************** !> ! Alternate version of [[feldg]] to be used with cartesian coordinates - subroutine feldc(me,v,b) + subroutine feldc(me, v, b) - class(shellig_type),intent(inout) :: me - real(wp),dimension(3),intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) + class(shellig_type), intent(inout) :: me + real(wp), dimension(3), intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) !! x-axis pointing to equator at 0 longitude !! y-axis pointing to equator at 90 long. !! z-axis pointing to north pole - real(wp),intent(out) :: b(3) !! field components - - real(wp) :: f , rq , s , t , x , xxx , y , yyy , z , zzz - integer :: i , ih , ihmax , il , imax , k , last , m - - xxx=v(1) - yyy=v(2) - zzz=v(3) - - rq=1.0_wp/(xxx*xxx+yyy*yyy+zzz*zzz) - me%xi = [xxx,yyy,zzz] * rq - - ihmax=me%nmax*me%nmax+1 - last=ihmax+me%nmax+me%nmax - imax=me%nmax+me%nmax-1 - do i=ihmax,last - me%h(i)=me%g(i) - end do - do k=1,3,2 - i=imax - ih=ihmax - do - il=ih-i - f=2.0_wp/real(i-k+2, wp) - x=me%xi(1)*f - y=me%xi(2)*f - z=me%xi(3)*(f+f) - i=i-2 - if ((i-1)>=0) then - if ((i-1)>0) then - do m=3,i,2 - me%h(il+m+1)=me%g(il+m+1)+z*me%h(ih+m+1)+x*(me%h(ih+m+3)-& - me%h(ih+m-1))-y*(me%h(ih+m+2)+me%h(ih+m-2)) - me%h(il+m)=me%g(il+m)+z*me%h(ih+m)+x*(me%h(ih+m+2)-& - me%h(ih+m-2))+y*(me%h(ih+m+3)+me%h(ih+m-1)) - end do - end if - me%h(il+2)=me%g(il+2)+z*me%h(ih+2)+x*me%h(ih+4)-y*(me%h(ih+3)+me%h(ih)) - me%h(il+1)=me%g(il+1)+z*me%h(ih+1)+y*me%h(ih+4)+x*(me%h(ih+3)-me%h(ih)) - end if - me%h(il)=me%g(il)+z*me%h(ih)+2.0_wp*(x*me%h(ih+1)+y*me%h(ih+2)) - ih=il - if (i= 0) then + if ((i - 1) > 0) then + do m = 3, i, 2 + me%h(il + m + 1) = me%g(il + m + 1) + z * me%h(ih + m + 1) + x * (me%h(ih + m + 3) - & + me%h(ih + m - 1)) - y * (me%h(ih + m + 2) + me%h(ih + m - 2)) + me%h(il + m) = me%g(il + m) + z * me%h(ih + m) + x * (me%h(ih + m + 2) - & + me%h(ih + m - 2)) + y * (me%h(ih + m + 3) + me%h(ih + m - 1)) + end do + end if + me%h(il + 2) = me%g(il + 2) + z * me%h(ih + 2) + x * me%h(ih + 4) - y * (me%h(ih + 3) + me%h(ih)) + me%h(il + 1) = me%g(il + 1) + z * me%h(ih + 1) + y * me%h(ih + 4) + x * (me%h(ih + 3) - me%h(ih)) + end if + me%h(il) = me%g(il) + z * me%h(ih) + 2.0_wp * (x * me%h(ih + 1) + y * me%h(ih + 2)) + ih = il + if (i < k) exit + end do + end do - s=0.5_wp*me%h(1)+2.0_wp*(me%h(2)*me%xi(3)+me%h(3)*me%xi(1)+me%h(4)*me%xi(2)) - t=(rq+rq)*sqrt(rq) + s = 0.5_wp * me%h(1) + 2.0_wp * (me%h(2) * me%xi(3) + me%h(3) * me%xi(1) + me%h(4) * me%xi(2)) + t = (rq + rq) * sqrt(rq) - b(1)=t*(me%h(3)-s*xxx) - b(2)=t*(me%h(4)-s*yyy) - b(3)=t*(me%h(2)-s*zzz) + b(1) = t * (me%h(3) - s * xxx) + b(2) = t * (me%h(4) - s * yyy) + b(3) = t * (me%h(2) - s * zzz) - end subroutine feldc + end subroutine feldc !***************************************************************************************** !> ! Used for `l` computation. - subroutine feldi(me) - - class(shellig_type),intent(inout) :: me - - real(wp) :: f , x , y , z - integer :: i , ih , ihmax , il , imax , k , last , m - - ihmax=me%nmax*me%nmax+1 - last=ihmax+me%nmax+me%nmax - imax=me%nmax+me%nmax-1 - do i=ihmax,last - me%h(i)=me%g(i) - end do - do k=1,3,2 - i=imax - ih=ihmax - do - il=ih-i - f=2.0_wp/real(i-k+2, wp) - x=me%xi(1)*f - y=me%xi(2)*f - z=me%xi(3)*(f+f) - i=i-2 - if ((i-1)>=0) then - if ((i-1)>0) then - do m=3,i,2 - me%h(il+m+1)=me%g(il+m+1)+z*me%h(ih+m+1)+x*(me%h(ih+m+3)-& - me%h(ih+m-1))-y*(me%h(ih+m+2)+me%h(ih+m-2)) - me%h(il+m)=me%g(il+m)+z*me%h(ih+m)+x*(me%h(ih+m+2)-& - me%h(ih+m-2))+y*(me%h(ih+m+3)+me%h(ih+m-1)) - end do - end if - me%h(il+2)=me%g(il+2)+z*me%h(ih+2)+x*me%h(ih+4)-y*(me%h(ih+3)+me%h(ih)) - me%h(il+1)=me%g(il+1)+z*me%h(ih+1)+y*me%h(ih+4)+x*(me%h(ih+3)-me%h(ih)) - end if - me%h(il)=me%g(il)+z*me%h(ih)+2.0_wp*(x*me%h(ih+1)+y*me%h(ih+2)) - ih=il - if (i= 0) then + if ((i - 1) > 0) then + do m = 3, i, 2 + me%h(il + m + 1) = me%g(il + m + 1) + z * me%h(ih + m + 1) + x * (me%h(ih + m + 3) - & + me%h(ih + m - 1)) - y * (me%h(ih + m + 2) + me%h(ih + m - 2)) + me%h(il + m) = me%g(il + m) + z * me%h(ih + m) + x * (me%h(ih + m + 2) - & + me%h(ih + m - 2)) + y * (me%h(ih + m + 3) + me%h(ih + m - 1)) + end do + end if + me%h(il + 2) = me%g(il + 2) + z * me%h(ih + 2) + x * me%h(ih + 4) - y * (me%h(ih + 3) + me%h(ih)) + me%h(il + 1) = me%g(il + 1) + z * me%h(ih + 1) + y * me%h(ih + 4) + x * (me%h(ih + 3) - me%h(ih)) + end if + me%h(il) = me%g(il) + z * me%h(ih) + 2.0_wp * (x * me%h(ih + 1) + y * me%h(ih + 2)) + ih = il + if (i < k) exit + end do + end do - end subroutine feldi + end subroutine feldi !***************************************************************************************** !> @@ -844,100 +844,100 @@ end subroutine feldi ! * updated to IGRF-2000 version -dkb- 5/31/2000 ! * updated to IGRF-2005 version -dkb- 3/24/2000 - subroutine feldcof(me,year,dimo) + subroutine feldcof(me, year, dimo) - class(shellig_type),intent(inout) :: me - real(wp),intent(in) :: year !! decimal year for which geomagnetic field is to + class(shellig_type), intent(inout) :: me + real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) - real(wp),intent(out) :: dimo !! geomagnetic dipol moment in gauss (normalized + real(wp), intent(out) :: dimo !! geomagnetic dipol moment in gauss (normalized !! to earth's radius) at the time (year) - real(wp) :: dte1 , dte2 , erad , gha(144) , sqrt2 - integer :: i , ier , j , l , m , n , iyea - character(len=:),allocatable :: fil2 - real(wp) :: x , f0 , f !! these were double precision in original + real(wp) :: dte1, dte2, erad, gha(144), sqrt2 + integer :: i, ier, j, l, m, n, iyea + character(len=:), allocatable :: fil2 + real(wp) :: x, f0, f !! these were double precision in original !! code while everything else was single precision - ! changed to conform with IGRF 45-95, also FILMOD, DTEMOD arrays +1 - character(len=filename_len),dimension(17),parameter :: filmod = [& - 'dgrf1945.dat ' , 'dgrf1950.dat ' , 'dgrf1955.dat ' , 'dgrf1960.dat ' , & - 'dgrf1965.dat ' , 'dgrf1970.dat ' , 'dgrf1975.dat ' , 'dgrf1980.dat ' , & - 'dgrf1985.dat ' , 'dgrf1990.dat ' , 'dgrf1995.dat ' , 'dgrf2000.dat ' , & - 'dgrf2005.dat ' , 'dgrf2010.dat ' , 'dgrf2015.dat ' , 'igrf2020.dat ' , & - 'igrf2020s.dat'] - real(wp),dimension(17),parameter :: dtemod = [1945.0_wp , 1950.0_wp , 1955.0_wp , & - 1960.0_wp , 1965.0_wp , 1970.0_wp , & - 1975.0_wp , 1980.0_wp , 1985.0_wp , & - 1990.0_wp , 1995.0_wp , 2000.0_wp , & - 2005.0_wp , 2010.0_wp , 2015.0_wp , & - 2020.0_wp , 2025.0_wp] - integer,parameter :: numye = size(dtemod)-1 ! number of 5-year priods represented by IGRF - integer,parameter :: is = 0 !! * is=0 for schmidt normalization + ! changed to conform with IGRF 45-95, also FILMOD, DTEMOD arrays +1 + character(len=filename_len), dimension(17), parameter :: filmod = [ & + 'dgrf1945.dat ', 'dgrf1950.dat ', 'dgrf1955.dat ', 'dgrf1960.dat ', & + 'dgrf1965.dat ', 'dgrf1970.dat ', 'dgrf1975.dat ', 'dgrf1980.dat ', & + 'dgrf1985.dat ', 'dgrf1990.dat ', 'dgrf1995.dat ', 'dgrf2000.dat ', & + 'dgrf2005.dat ', 'dgrf2010.dat ', 'dgrf2015.dat ', 'igrf2020.dat ', & + 'igrf2020s.dat'] + real(wp), dimension(17), parameter :: dtemod = [1945.0_wp, 1950.0_wp, 1955.0_wp, & + 1960.0_wp, 1965.0_wp, 1970.0_wp, & + 1975.0_wp, 1980.0_wp, 1985.0_wp, & + 1990.0_wp, 1995.0_wp, 2000.0_wp, & + 2005.0_wp, 2010.0_wp, 2015.0_wp, & + 2020.0_wp, 2025.0_wp] + integer, parameter :: numye = size(dtemod) - 1 ! number of 5-year priods represented by IGRF + integer, parameter :: is = 0 !! * is=0 for schmidt normalization !! * is=1 gauss normalization - logical :: read_file - - !-- determine igrf-years for input-year - me%time = year - iyea = int(year/5.0_wp)*5 - read_file = iyea /= me%iyea ! if we have to read the file - me%iyea = iyea - l = (me%iyea-1945)/5 + 1 - if ( l<1 ) l = 1 - if ( l>numye ) l = numye - dte1 = dtemod(l) - me%name = me%get_data_file_dir() // trim(filmod(l)) - dte2 = dtemod(l+1) - fil2 = me%get_data_file_dir() // trim(filmod(l+1)) - if (read_file) then - ! get igrf coefficients for the boundary years - ! [if they have not ready been loaded] - call getshc(me%name,me%nmax1,erad,me%g,ier) - if ( ier/=0 ) error stop 'error reading file: '//trim(me%name) - me%g_cache = me%g ! because it is modified below, we have to cache the original values from the file - call getshc(fil2,me%nmax2,erad,me%gh2,ier) - if ( ier/=0 ) error stop 'error reading file: '//trim(fil2) - else - me%g = me%g_cache - end if - !-- determine igrf coefficients for year - if ( l<=numye-1 ) then - call intershc(year,dte1,me%nmax1,me%g,dte2,me%nmax2,me%gh2,me%nmax,gha) - else - call extrashc(year,dte1,me%nmax1,me%g,me%nmax2,me%gh2,me%nmax,gha) - endif - !-- determine magnetic dipol moment and coeffiecients g - f0 = 0.0_wp - do j = 1 , 3 - f = gha(j)*1.0e-5_wp - f0 = f0 + f*f - enddo - dimo = sqrt(f0) - - me%g(1) = 0.0_wp - i = 2 - f0 = 1.0e-5_wp - if ( is==0 ) f0 = -f0 - sqrt2 = sqrt(2.0_wp) - - do n = 1 , me%nmax - x = n - f0 = f0*x*x/(4.0_wp*x-2.0_wp) - if ( is==0 ) f0 = f0*(2.0_wp*x-1.0_wp)/x - f = f0*0.5_wp - if ( is==0 ) f = f*sqrt2 - me%g(i) = gha(i-1)*f0 - i = i + 1 - do m = 1 , n - f = f*(x+m)/(x-m+1.0_wp) - if ( is==0 ) f = f*sqrt((x-m+1.0_wp)/(x+m)) - me%g(i) = gha(i-1)*f - me%g(i+1) = gha(i)*f - i = i + 2 - enddo - enddo - -end subroutine feldcof + logical :: read_file + + !-- determine igrf-years for input-year + me%time = year + iyea = int(year / 5.0_wp) * 5 + read_file = iyea /= me%iyea ! if we have to read the file + me%iyea = iyea + l = (me%iyea - 1945) / 5 + 1 + if (l < 1) l = 1 + if (l > numye) l = numye + dte1 = dtemod(l) + me%name = me%get_data_file_dir()//trim(filmod(l)) + dte2 = dtemod(l + 1) + fil2 = me%get_data_file_dir()//trim(filmod(l + 1)) + if (read_file) then + ! get igrf coefficients for the boundary years + ! [if they have not ready been loaded] + call getshc(me%name, me%nmax1, erad, me%g, ier) + if (ier /= 0) error stop 'error reading file: '//trim(me%name) + me%g_cache = me%g ! because it is modified below, we have to cache the original values from the file + call getshc(fil2, me%nmax2, erad, me%gh2, ier) + if (ier /= 0) error stop 'error reading file: '//trim(fil2) + else + me%g = me%g_cache + end if + !-- determine igrf coefficients for year + if (l <= numye - 1) then + call intershc(year, dte1, me%nmax1, me%g, dte2, me%nmax2, me%gh2, me%nmax, gha) + else + call extrashc(year, dte1, me%nmax1, me%g, me%nmax2, me%gh2, me%nmax, gha) + end if + !-- determine magnetic dipol moment and coeffiecients g + f0 = 0.0_wp + do j = 1, 3 + f = gha(j) * 1.0e-5_wp + f0 = f0 + f * f + end do + dimo = sqrt(f0) + + me%g(1) = 0.0_wp + i = 2 + f0 = 1.0e-5_wp + if (is == 0) f0 = -f0 + sqrt2 = sqrt(2.0_wp) + + do n = 1, me%nmax + x = n + f0 = f0 * x * x / (4.0_wp * x - 2.0_wp) + if (is == 0) f0 = f0 * (2.0_wp * x - 1.0_wp) / x + f = f0 * 0.5_wp + if (is == 0) f = f * sqrt2 + me%g(i) = gha(i - 1) * f0 + i = i + 1 + do m = 1, n + f = f * (x + m) / (x - m + 1.0_wp) + if (is == 0) f = f * sqrt((x - m + 1.0_wp) / (x + m)) + me%g(i) = gha(i - 1) * f + me%g(i + 1) = gha(i) * f + i = i + 2 + end do + end do + + end subroutine feldcof !***************************************************************************************** !> @@ -948,82 +948,82 @@ end subroutine feldcof ! * Version 1.01, A. Zunde, USGS, MS 964, ! Box 25046 Federal Center, Denver, CO 80225 -subroutine getshc(Fspec,Nmax,Erad,Gh,Ier) + subroutine getshc(Fspec, Nmax, Erad, Gh, Ier) - character(len=*),intent(in) :: Fspec !! File specification - integer,intent(out) :: Nmax !! Maximum degree and order of model - real(wp),intent(out) :: Erad !! Earth's radius associated with the spherical + character(len=*), intent(in) :: Fspec !! File specification + integer, intent(out) :: Nmax !! Maximum degree and order of model + real(wp), intent(out) :: Erad !! Earth's radius associated with the spherical !! harmonic coefficients, in the same units as !! elevation - real(wp),dimension(*),intent(out) :: Gh !! Schmidt quasi-normal internal spherical + real(wp), dimension(*), intent(out) :: Gh !! Schmidt quasi-normal internal spherical !! harmonic coefficients - integer,intent(out) :: Ier !! Error number: + integer, intent(out) :: Ier !! Error number: !! !! * 0, no error !! * -2, records out of order !! * FORTRAN run-time error number - integer :: iu !! logical unit number - real(wp) :: g , h - integer :: i , m , mm , n , nn - - read_file : block - ! --------------------------------------------------------------- - ! Open coefficient file. Read past first header record. - ! Read degree and order of model and Earth's radius. - ! --------------------------------------------------------------- - OPEN (newunit=Iu,FILE=Fspec,STATUS='OLD',IOSTAT=Ier) - if (Ier/=0) then - write(*,*) 'Error opening file: '//trim(fspec) - exit read_file - end if - READ (Iu,*,IOSTAT=Ier) - if (Ier/=0) exit read_file - READ (Iu,*,IOSTAT=Ier) Nmax , Erad - if (Ier/=0) exit read_file - - ! --------------------------------------------------------------- - ! Read the coefficient file, arranged as follows: - ! - ! N M G H - ! ---------------------- - ! / 1 0 GH(1) - - ! / 1 1 GH(2) GH(3) - ! / 2 0 GH(4) - - ! / 2 1 GH(5) GH(6) - ! NMAX*(NMAX+3)/2 / 2 2 GH(7) GH(8) - ! records \ 3 0 GH(9) - - ! \ . . . . - ! \ . . . . - ! NMAX*(NMAX+2) \ . . . . - ! elements in GH \ NMAX NMAX . . - ! - ! N and M are, respectively, the degree and order of the - ! coefficient. - ! --------------------------------------------------------------- - i = 0 - main: DO nn = 1 , Nmax - DO mm = 0 , nn - READ (Iu,*,IOSTAT=Ier) n , m , g , h - if (Ier/=0) exit main - IF ( nn/=n .OR. mm/=m ) THEN - Ier = -2 - EXIT main - ENDIF - i = i + 1 - Gh(i) = g - IF ( m/=0 ) THEN - i = i + 1 - Gh(i) = h - ENDIF - ENDDO - ENDDO main + integer :: iu !! logical unit number + real(wp) :: g, h + integer :: i, m, mm, n, nn + + read_file: block + ! --------------------------------------------------------------- + ! Open coefficient file. Read past first header record. + ! Read degree and order of model and Earth's radius. + ! --------------------------------------------------------------- + open (newunit=Iu, FILE=Fspec, STATUS='OLD', IOSTAT=Ier) + if (Ier /= 0) then + write (*, *) 'Error opening file: '//trim(fspec) + exit read_file + end if + read (Iu, *, IOSTAT=Ier) + if (Ier /= 0) exit read_file + read (Iu, *, IOSTAT=Ier) Nmax, Erad + if (Ier /= 0) exit read_file + + ! --------------------------------------------------------------- + ! Read the coefficient file, arranged as follows: + ! + ! N M G H + ! ---------------------- + ! / 1 0 GH(1) - + ! / 1 1 GH(2) GH(3) + ! / 2 0 GH(4) - + ! / 2 1 GH(5) GH(6) + ! NMAX*(NMAX+3)/2 / 2 2 GH(7) GH(8) + ! records \ 3 0 GH(9) - + ! \ . . . . + ! \ . . . . + ! NMAX*(NMAX+2) \ . . . . + ! elements in GH \ NMAX NMAX . . + ! + ! N and M are, respectively, the degree and order of the + ! coefficient. + ! --------------------------------------------------------------- + i = 0 + main: do nn = 1, Nmax + do mm = 0, nn + read (Iu, *, IOSTAT=Ier) n, m, g, h + if (Ier /= 0) exit main + if (nn /= n .or. mm /= m) then + Ier = -2 + exit main + end if + i = i + 1 + Gh(i) = g + if (m /= 0) then + i = i + 1 + Gh(i) = h + end if + end do + end do main - end block read_file + end block read_file - CLOSE (Iu) + close (Iu) -END subroutine getshc + end subroutine getshc !***************************************************************************************** !> @@ -1041,47 +1041,47 @@ END subroutine getshc ! * Version 1.01, A. Zunde ! USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 -subroutine intershc(date,dte1,nmax1,gh1,dte2,nmax2,gh2,nmax,gh) - - real(wp),intent(in) :: date !! Date of resulting model (in decimal year) - real(wp),intent(in) :: dte1 !! Date of earlier model - integer,intent(in) :: nmax1 !! Maximum degree and order of earlier model - real(wp),intent(in) :: gh1(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of earlier model - real(wp),intent(in) :: dte2 !! Date of later model - integer,intent(in) :: nmax2 !! Maximum degree and order of later model - real(wp),intent(in) :: gh2(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of later model - real(wp),intent(out) :: gh(*) !! Coefficients of resulting model - integer,intent(out) :: nmax !! Maximum degree and order of resulting model - - real(wp) :: factor - integer :: i , k , l - - factor = (date-dte1)/(dte2-dte1) - - if ( nmax1==nmax2 ) then - k = nmax1*(nmax1+2) - nmax = nmax1 - elseif ( nmax1>nmax2 ) then - k = nmax2*(nmax2+2) - l = nmax1*(nmax1+2) - do i = k + 1 , l - gh(i) = gh1(i) + factor*(-gh1(i)) - enddo - nmax = nmax1 - else - k = nmax1*(nmax1+2) - l = nmax2*(nmax2+2) - do i = k + 1 , l - gh(i) = factor*gh2(i) - enddo - nmax = nmax2 - endif - - do i = 1 , k - gh(i) = gh1(i) + factor*(gh2(i)-gh1(i)) - enddo - -end subroutine intershc + subroutine intershc(date, dte1, nmax1, gh1, dte2, nmax2, gh2, nmax, gh) + + real(wp), intent(in) :: date !! Date of resulting model (in decimal year) + real(wp), intent(in) :: dte1 !! Date of earlier model + integer, intent(in) :: nmax1 !! Maximum degree and order of earlier model + real(wp), intent(in) :: gh1(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of earlier model + real(wp), intent(in) :: dte2 !! Date of later model + integer, intent(in) :: nmax2 !! Maximum degree and order of later model + real(wp), intent(in) :: gh2(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of later model + real(wp), intent(out) :: gh(*) !! Coefficients of resulting model + integer, intent(out) :: nmax !! Maximum degree and order of resulting model + + real(wp) :: factor + integer :: i, k, l + + factor = (date - dte1) / (dte2 - dte1) + + if (nmax1 == nmax2) then + k = nmax1 * (nmax1 + 2) + nmax = nmax1 + elseif (nmax1 > nmax2) then + k = nmax2 * (nmax2 + 2) + l = nmax1 * (nmax1 + 2) + do i = k + 1, l + gh(i) = gh1(i) + factor * (-gh1(i)) + end do + nmax = nmax1 + else + k = nmax1 * (nmax1 + 2) + l = nmax2 * (nmax2 + 2) + do i = k + 1, l + gh(i) = factor * gh2(i) + end do + nmax = nmax2 + end if + + do i = 1, k + gh(i) = gh1(i) + factor * (gh2(i) - gh1(i)) + end do + + end subroutine intershc !***************************************************************************************** !> @@ -1099,77 +1099,77 @@ end subroutine intershc ! * Version 1.01, A. Zunde ! USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 -subroutine extrashc(date,dte1,nmax1,gh1,nmax2,gh2,nmax,gh) - - real(wp),intent(in) :: date !! Date of resulting model (in decimal year) - real(wp),intent(in) :: dte1 !! Date of base model - integer,intent(in) :: nmax1 !! Maximum degree and order of base model - real(wp),intent(in) :: gh1(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of base model - integer,intent(in) :: nmax2 !! Maximum degree and order of rate-of-change model - real(wp),intent(in) :: gh2(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of rate-of-change model - real(wp),intent(out) :: gh(*) !! Coefficients of resulting model - integer,intent(out) :: nmax !! Maximum degree and order of resulting model - - real(wp) :: factor - integer :: i , k , l - - factor = (date-dte1) - - if ( nmax1==nmax2 ) then - k = nmax1*(nmax1+2) - nmax = nmax1 - elseif ( nmax1>nmax2 ) then - k = nmax2*(nmax2+2) - l = nmax1*(nmax1+2) - do i = k + 1 , l - gh(i) = gh1(i) - enddo - nmax = nmax1 - else - k = nmax1*(nmax1+2) - l = nmax2*(nmax2+2) - do i = k + 1 , l - gh(i) = factor*gh2(i) - enddo - nmax = nmax2 - endif - - do i = 1 , k - gh(i) = gh1(i) + factor*gh2(i) - enddo - -end subroutine extrashc + subroutine extrashc(date, dte1, nmax1, gh1, nmax2, gh2, nmax, gh) + + real(wp), intent(in) :: date !! Date of resulting model (in decimal year) + real(wp), intent(in) :: dte1 !! Date of base model + integer, intent(in) :: nmax1 !! Maximum degree and order of base model + real(wp), intent(in) :: gh1(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of base model + integer, intent(in) :: nmax2 !! Maximum degree and order of rate-of-change model + real(wp), intent(in) :: gh2(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of rate-of-change model + real(wp), intent(out) :: gh(*) !! Coefficients of resulting model + integer, intent(out) :: nmax !! Maximum degree and order of resulting model + + real(wp) :: factor + integer :: i, k, l + + factor = (date - dte1) + + if (nmax1 == nmax2) then + k = nmax1 * (nmax1 + 2) + nmax = nmax1 + elseif (nmax1 > nmax2) then + k = nmax2 * (nmax2 + 2) + l = nmax1 * (nmax1 + 2) + do i = k + 1, l + gh(i) = gh1(i) + end do + nmax = nmax1 + else + k = nmax1 * (nmax1 + 2) + l = nmax2 * (nmax2 + 2) + do i = k + 1, l + gh(i) = factor * gh2(i) + end do + nmax = nmax2 + end if + + do i = 1, k + gh(i) = gh1(i) + factor * gh2(i) + end do + + end subroutine extrashc !***************************************************************************************** !> ! geodetic to scaled cartesian coordinates -pure function geo_to_cart(glat,glon,alt) result(x) + pure function geo_to_cart(glat, glon, alt) result(x) - real(wp),intent(in) :: glat !! geodetic latitude in degrees (north) - real(wp),intent(in) :: glon !! geodetic longitude in degrees (east) - real(wp),intent(in) :: alt !! altitude in km above sea level - real(wp),dimension(3) :: x !! cartesian coordinates in earth radii (6371.2 km) + real(wp), intent(in) :: glat !! geodetic latitude in degrees (north) + real(wp), intent(in) :: glon !! geodetic longitude in degrees (east) + real(wp), intent(in) :: alt !! altitude in km above sea level + real(wp), dimension(3) :: x !! cartesian coordinates in earth radii (6371.2 km) !! !! * x-axis pointing to equator at 0 longitude !! * y-axis pointing to equator at 90 long. !! * z-axis pointing to north pole - real(wp) :: rlat !! latitude in radians - real(wp) :: rlon !! longitude in radians - real(wp) :: d, rho + real(wp) :: rlat !! latitude in radians + real(wp) :: rlon !! longitude in radians + real(wp) :: d, rho - ! deg to radians: - rlat = glat*umr - rlon = glon*umr + ! deg to radians: + rlat = glat * umr + rlon = glon * umr - ! JW : it's weird that ct is sin, and st is cos...it was like that in the original code - associate (ct => sin(rlat), st => cos(rlat), cp => cos(rlon), sp => sin(rlon)) - d = sqrt(aquad-(aquad-bquad)*ct*ct) - rho = (alt+aquad/d)*st/era - x = [rho*cp, rho*sp, (alt+bquad/d)*ct/era] - end associate + ! JW : it's weird that ct is sin, and st is cos...it was like that in the original code + associate (ct => sin(rlat), st => cos(rlat), cp => cos(rlon), sp => sin(rlon)) + d = sqrt(aquad - (aquad - bquad) * ct * ct) + rho = (alt + aquad / d) * st / era + x = [rho * cp, rho * sp, (alt + bquad / d) * ct / era] + end associate -end function geo_to_cart + end function geo_to_cart -end module SHELLIG_module \ No newline at end of file +end module SHELLIG_module diff --git a/src/trmfun.f90 b/src/trmfun.f90 index feacd6d..19dfac3 100644 --- a/src/trmfun.f90 +++ b/src/trmfun.f90 @@ -7,121 +7,121 @@ module trmfun_module - use radbelt_kinds_module + use radbelt_kinds_module - implicit none + implicit none - private + private - character(len=10),dimension(4),parameter :: mname = [ 'ae8min.asc' , & - 'ae8max.asc' , & - 'ap8min.asc' , & - 'ap8max.asc'] !! data files available + character(len=10), dimension(4), parameter :: mname = ['ae8min.asc', & + 'ae8max.asc', & + 'ap8min.asc', & + 'ap8max.asc'] !! data files available - type,public :: trm_type + type, public :: trm_type !! main class for the `aep8` model - private + private - character(len=:),allocatable :: aep8_dir !! directory containing the data files + character(len=:), allocatable :: aep8_dir !! directory containing the data files - ! data read from the files: - character(len=:),allocatable :: file_loaded !! the file that has been loaded - integer,dimension(8) :: ihead = 0 - integer,dimension(:),allocatable :: map + ! data read from the files: + character(len=:), allocatable :: file_loaded !! the file that has been loaded + integer, dimension(8) :: ihead = 0 + integer, dimension(:), allocatable :: map - real(wp) :: fistep = 0.0_wp !! the stepsize for the parameterization of the logarithm of flux. + real(wp) :: fistep = 0.0_wp !! the stepsize for the parameterization of the logarithm of flux. !! formerly stored in common block `tra2` - ! formerly saved variables in trara1: - real(wp) :: f1 = 1.001_wp - real(wp) :: f2 = 1.002_wp + ! formerly saved variables in trara1: + real(wp) :: f1 = 1.001_wp + real(wp) :: f2 = 1.002_wp - contains - private - procedure,public :: aep8 !! main routine - procedure,public :: trara1, trara2 !! low-level routine - procedure,public :: set_data_file_dir, get_data_file_dir - end type trm_type + contains + private + procedure, public :: aep8 !! main routine + procedure, public :: trara1, trara2 !! low-level routine + procedure, public :: set_data_file_dir, get_data_file_dir + end type trm_type - contains +contains !***************************************************************************************** !> ! Set the directory containing the data files. - subroutine set_data_file_dir(me,dir) - class(trm_type),intent(inout) :: me - character(len=*),intent(in) :: dir - me%aep8_dir = trim(dir) - end subroutine set_data_file_dir + subroutine set_data_file_dir(me, dir) + class(trm_type), intent(inout) :: me + character(len=*), intent(in) :: dir + me%aep8_dir = trim(dir) + end subroutine set_data_file_dir !***************************************************************************************** !> ! Get the directory containing the data files. - function get_data_file_dir(me) result(dir) - class(trm_type),intent(in) :: me - character(len=:),allocatable :: dir - if (allocated(me%aep8_dir)) then - dir = trim(me%aep8_dir) // '/' - else - dir = 'data/aep8/' ! default - end if - end function get_data_file_dir + function get_data_file_dir(me) result(dir) + class(trm_type), intent(in) :: me + character(len=:), allocatable :: dir + if (allocated(me%aep8_dir)) then + dir = trim(me%aep8_dir)//'/' + else + dir = 'data/aep8/' ! default + end if + end function get_data_file_dir !***************************************************************************************** !> ! Main wrapper for the radiation model. ! Reads the coefficient file and calls the low-level routine. - subroutine aep8(me,e,l,bb0,imname,flux) - - class(trm_type),intent(inout) :: me - - real(wp),intent(in) :: e - real(wp),intent(in) :: l - real(wp),intent(in) :: bb0 - integer,intent(in) :: imname !! which model to load (index in `mname` array) - real(wp),intent(out) :: flux - - real(wp) :: ee(1), f(1) !! temp variables - integer :: i , ierr, iuaeap , nmap - character(len=:),allocatable :: name - logical :: load_file - - name = me%get_data_file_dir() // trim(mname(Imname)) ! the file to load - - ! JW : do we need to reset some or all of these ? - me%fistep = 0.0_wp - me%f1 = 1.001_wp - me%f2 = 1.002_wp - - ! check to see if this file has already been loaded - ! [the class can store one file at a time] - load_file = .true. - if (allocated(me%file_loaded)) then - if (name == me%file_loaded) load_file = .false. - end if - - if (load_file) then - open (newunit = iuaeap,file=name,status='OLD',iostat=ierr,form='FORMATTED') - if ( ierr/=0 ) then - error stop 'error reading '//trim(name) - end if - read (iuaeap,'(1X,12I6)') me%ihead - nmap = me%ihead(8) - allocate(me%map(nmap)) - read (iuaeap,'(1X,12I6)') (me%map(i),i=1,nmap) - close (iuaeap) - me%file_loaded = trim(name) - end if - - ee(1) = e - call me%trara1(me%ihead,me%map,L,Bb0,ee,f,1) - flux = f(1) - IF ( Flux>0.0_wp ) Flux = 10.0_wp**Flux - - end subroutine aep8 + subroutine aep8(me, e, l, bb0, imname, flux) + + class(trm_type), intent(inout) :: me + + real(wp), intent(in) :: e + real(wp), intent(in) :: l + real(wp), intent(in) :: bb0 + integer, intent(in) :: imname !! which model to load (index in `mname` array) + real(wp), intent(out) :: flux + + real(wp) :: ee(1), f(1) !! temp variables + integer :: i, ierr, iuaeap, nmap + character(len=:), allocatable :: name + logical :: load_file + + name = me%get_data_file_dir()//trim(mname(Imname)) ! the file to load + + ! JW : do we need to reset some or all of these ? + me%fistep = 0.0_wp + me%f1 = 1.001_wp + me%f2 = 1.002_wp + + ! check to see if this file has already been loaded + ! [the class can store one file at a time] + load_file = .true. + if (allocated(me%file_loaded)) then + if (name == me%file_loaded) load_file = .false. + end if + + if (load_file) then + open (newunit=iuaeap, file=name, status='OLD', iostat=ierr, form='FORMATTED') + if (ierr /= 0) then + error stop 'error reading '//trim(name) + end if + read (iuaeap, '(1X,12I6)') me%ihead + nmap = me%ihead(8) + allocate (me%map(nmap)) + read (iuaeap, '(1X,12I6)') (me%map(i), i=1, nmap) + close (iuaeap) + me%file_loaded = trim(name) + end if + + ee(1) = e + call me%trara1(me%ihead, me%map, L, Bb0, ee, f, 1) + flux = f(1) + if (Flux > 0.0_wp) Flux = 10.0_wp**Flux + + end subroutine aep8 !***************************************************************************************** !***************************************************************************************** @@ -130,114 +130,114 @@ end subroutine aep8 ! strength and l-value. function [[trara2]] is used to interpolate in ! b-l-space. - subroutine trara1(me,descr,map,fl,bb0,e,f,n) + subroutine trara1(me, descr, map, fl, bb0, e, f, n) - class(trm_type),intent(inout) :: me - integer,intent(in) :: n !! number of energies - integer,intent(in) :: descr(8) !! header of specified trapped radition model - real(wp),intent(in) :: e(n) !! array of energies in mev - real(wp),intent(in) :: fl !! l-value - real(wp),intent(in) :: bb0 !! =b/b0 magnetic field strength normalized + class(trm_type), intent(inout) :: me + integer, intent(in) :: n !! number of energies + integer, intent(in) :: descr(8) !! header of specified trapped radition model + real(wp), intent(in) :: e(n) !! array of energies in mev + real(wp), intent(in) :: fl !! l-value + real(wp), intent(in) :: bb0 !! =b/b0 magnetic field strength normalized !! to field strength at magnetic equator - integer,intent(in) :: map(*) !! map of trapped radition model + integer, intent(in) :: map(*) !! map of trapped radition model !! (descr and map are explained at the begin !! of the main program model) - real(wp),intent(out) :: f(n) !! decadic logarithm of integral fluxes in + real(wp), intent(out) :: f(n) !! decadic logarithm of integral fluxes in !! particles/(cm*cm*sec) - real(wp) :: e0 , e1 , e2 , escale , f0 , fscale , xnl - real(wp) :: bb0_ !! local copy of `bb0`. in the original code + real(wp) :: e0, e1, e2, escale, f0, fscale, xnl + real(wp) :: bb0_ !! local copy of `bb0`. in the original code !! this was modified by this routine. !! added this so `bb0` could be `intent(in)` - integer :: i0 , i1 , i2 , i3 , ie , l3 , nb , nl - logical :: s0 , s1 , s2 - - e0 = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings - f0 = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings - i0 = 0 ! to avoid -Wmaybe-uninitialized warnings - s0 = .false. ! to avoid -Wmaybe-uninitialized warnings -- but not sure what default value here should be ! -JW - - bb0_ = bb0 - me%fistep = descr(7)/descr(2) - escale = descr(4) - fscale = descr(7) - xnl = min(15.6_wp,abs(fl)) - nl = int(xnl*descr(5)) - if ( bb0_<1.0_wp ) bb0_ = 1.0_wp - nb = int((bb0_-1.0_wp)*descr(6)) - - ! i2 is the number of elements in the flux map for the first energy. - ! i3 is the index of the last element of the second energy map. - ! l3 is the length of the map for the third energy. - ! e1 is the energy of the first energy map (unscaled) - ! e2 is the energy of the second energy map (unscaled) - i1 = 0 - i2 = map(1) - i3 = i2 + map(i2+1) - l3 = map(i3+1) - e1 = map(i1+2)/escale - e2 = map(i2+2)/escale - - ! s0, s1, s2 are logical variables which indicate whether the flux for - ! a particular e, b, l point has already been found in a previous call - ! to function trara2. if not, s.. =.true. - s1 = .true. - s2 = .true. - - ! energy loop - - do ie = 1 , n - - ! for each energy e(i) find the successive energies e0,e1,e2 in - ! model map, which obey e0 < e1 < e(i) < e2 . - - do while ( (e(ie)>e2) .and. (l3/=0) ) - i0 = i1 - i1 = i2 - i2 = i3 - i3 = i3 + l3 - l3 = map(i3+1) - e0 = e1 - e1 = e2 - e2 = map(i2+2)/escale - s0 = s1 - s1 = s2 - s2 = .true. - f0 = me%f1 - me%f1 = me%f2 - enddo - - ! call trara2 to interpolate the flux-maps for e1,e2 in l-b/b0- - ! space to find fluxes f1,f2 [if they have not already been - ! calculated for a previous e(i)]. - - if ( s1 ) me%f1 = me%trara2(map(i1+3),nl,nb)/fscale - if ( s2 ) me%f2 = me%trara2(map(i2+3),nl,nb)/fscale - s1 = .false. - s2 = .false. - - ! finally, interpolate in energy. - - f(ie) = me%f1 + (me%f2-me%f1)*(e(ie)-e1)/(e2-e1) - if ( me%f2<=0.0_wp ) then - if ( i1/=0 ) then - ! --------- special interpolation --------------------------------- - ! if the flux for the second energy cannot be found (i.e. f2=0.0), - ! and the zeroth energy map has been defined (i.e. i1 not equal 0), - ! then interpolate using the flux maps for the zeroth and first - ! energy and choose the minimum of this interpolations and the - ! interpolation that was done with f2=0. - if ( s0 ) f0 = me%trara2(map(i0+3),nl,nb)/fscale - s0 = .false. - f(ie) = min(f(ie),f0+(me%f1-f0)*(e(ie)-e0)/(e1-e0)) - endif - endif - - ! the logarithmic flux is always kept greater or equal zero. - - f(ie) = max(f(ie),0.0_wp) - enddo -end subroutine trara1 + integer :: i0, i1, i2, i3, ie, l3, nb, nl + logical :: s0, s1, s2 + + e0 = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings + f0 = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings + i0 = 0 ! to avoid -Wmaybe-uninitialized warnings + s0 = .false. ! to avoid -Wmaybe-uninitialized warnings -- but not sure what default value here should be ! -JW + + bb0_ = bb0 + me%fistep = descr(7) / descr(2) + escale = descr(4) + fscale = descr(7) + xnl = min(15.6_wp, abs(fl)) + nl = int(xnl * descr(5)) + if (bb0_ < 1.0_wp) bb0_ = 1.0_wp + nb = int((bb0_ - 1.0_wp) * descr(6)) + + ! i2 is the number of elements in the flux map for the first energy. + ! i3 is the index of the last element of the second energy map. + ! l3 is the length of the map for the third energy. + ! e1 is the energy of the first energy map (unscaled) + ! e2 is the energy of the second energy map (unscaled) + i1 = 0 + i2 = map(1) + i3 = i2 + map(i2 + 1) + l3 = map(i3 + 1) + e1 = map(i1 + 2) / escale + e2 = map(i2 + 2) / escale + + ! s0, s1, s2 are logical variables which indicate whether the flux for + ! a particular e, b, l point has already been found in a previous call + ! to function trara2. if not, s.. =.true. + s1 = .true. + s2 = .true. + + ! energy loop + + do ie = 1, n + + ! for each energy e(i) find the successive energies e0,e1,e2 in + ! model map, which obey e0 < e1 < e(i) < e2 . + + do while ((e(ie) > e2) .and. (l3 /= 0)) + i0 = i1 + i1 = i2 + i2 = i3 + i3 = i3 + l3 + l3 = map(i3 + 1) + e0 = e1 + e1 = e2 + e2 = map(i2 + 2) / escale + s0 = s1 + s1 = s2 + s2 = .true. + f0 = me%f1 + me%f1 = me%f2 + end do + + ! call trara2 to interpolate the flux-maps for e1,e2 in l-b/b0- + ! space to find fluxes f1,f2 [if they have not already been + ! calculated for a previous e(i)]. + + if (s1) me%f1 = me%trara2(map(i1 + 3), nl, nb) / fscale + if (s2) me%f2 = me%trara2(map(i2 + 3), nl, nb) / fscale + s1 = .false. + s2 = .false. + + ! finally, interpolate in energy. + + f(ie) = me%f1 + (me%f2 - me%f1) * (e(ie) - e1) / (e2 - e1) + if (me%f2 <= 0.0_wp) then + if (i1 /= 0) then + ! --------- special interpolation --------------------------------- + ! if the flux for the second energy cannot be found (i.e. f2=0.0), + ! and the zeroth energy map has been defined (i.e. i1 not equal 0), + ! then interpolate using the flux maps for the zeroth and first + ! energy and choose the minimum of this interpolations and the + ! interpolation that was done with f2=0. + if (s0) f0 = me%trara2(map(i0 + 3), nl, nb) / fscale + s0 = .false. + f(ie) = min(f(ie), f0 + (me%f1 - f0) * (e(ie) - e0) / (e1 - e0)) + end if + end if + + ! the logarithmic flux is always kept greater or equal zero. + + f(ie) = max(f(ie), 0.0_wp) + end do + end subroutine trara1 !***************************************************************************************** !> @@ -248,258 +248,258 @@ end subroutine trara1 ! see main program 'model' for explanation of map format ! scaling factors. -function trara2(me,map,il,ib) + function trara2(me, map, il, ib) - class(trm_type),intent(inout) :: me - integer,intent(in) :: map(*) !! is sub-map (for specific energy) of + class(trm_type), intent(inout) :: me + integer, intent(in) :: map(*) !! is sub-map (for specific energy) of !! trapped radiation model map - integer,intent(in) :: il !! scaled l-value - integer,intent(in) :: ib !! scaled b/b0-1 - real(wp) :: trara2 !! scaled logarithm of particle flux - - real(wp) :: dfl , fincr1 , fincr2 , fistep , fkb , fkb1 , fkb2 , fkbj1 , fkbj2 , & - fkbm , fll1 , fll2 , flog , flog1 , flog2 , flogm , & - fnb , fnl , sl1 , sl2 - integer :: i1 , i2 , itime , j1 , j2 , kt , l1 , l2 - logical :: dummy - - fistep = me%fistep - - !******** - ! to avoid -Wmaybe-uninitialized warning - dfl = 0.0_wp - fincr1 = 0.0_wp - fincr2 = 0.0_wp - fkb = 0.0_wp - fkb1 = 0.0_wp - fkb2 = 0.0_wp - fkbm = 0.0_wp - flog = 0.0_wp - flog1 = 0.0_wp - flog2 = 0.0_wp - flogm = 0.0_wp - fnb = 0.0_wp - fnl = 0.0_wp - sl2 = 0.0_wp - i1 = 0 - i2 = 0 - itime = 0 - j2 = 0 - l1 = 0 - l2 = 0 - !******** - - ! these are recursive functions that - ! replace the gotos in the original code - call task1(dummy) - - contains - - recursive subroutine task1(done) - logical,intent(out) :: done - done = .false. - fnl = il - fnb = ib - itime = 0 - i2 = 0 - do - ! find consecutive sub-sub-maps for scaled l-values ls1,ls2, - ! with il less or equal ls2. l1,l2 are lengths of sub-sub-maps. - ! i1,i2 are indeces of first elements minus 1. - l2 = map(i2+1) - if ( map(i2+2)<=il ) then + integer, intent(in) :: il !! scaled l-value + integer, intent(in) :: ib !! scaled b/b0-1 + real(wp) :: trara2 !! scaled logarithm of particle flux + + real(wp) :: dfl, fincr1, fincr2, fistep, fkb, fkb1, fkb2, fkbj1, fkbj2, & + fkbm, fll1, fll2, flog, flog1, flog2, flogm, & + fnb, fnl, sl1, sl2 + integer :: i1, i2, itime, j1, j2, kt, l1, l2 + logical :: dummy + + fistep = me%fistep + + !******** + ! to avoid -Wmaybe-uninitialized warning + dfl = 0.0_wp + fincr1 = 0.0_wp + fincr2 = 0.0_wp + fkb = 0.0_wp + fkb1 = 0.0_wp + fkb2 = 0.0_wp + fkbm = 0.0_wp + flog = 0.0_wp + flog1 = 0.0_wp + flog2 = 0.0_wp + flogm = 0.0_wp + fnb = 0.0_wp + fnl = 0.0_wp + sl2 = 0.0_wp + i1 = 0 + i2 = 0 + itime = 0 + j2 = 0 + l1 = 0 + l2 = 0 + !******** + + ! these are recursive functions that + ! replace the gotos in the original code + call task1(dummy) + + contains + + recursive subroutine task1(done) + logical, intent(out) :: done + done = .false. + fnl = il + fnb = ib + itime = 0 + i2 = 0 + do + ! find consecutive sub-sub-maps for scaled l-values ls1,ls2, + ! with il less or equal ls2. l1,l2 are lengths of sub-sub-maps. + ! i1,i2 are indeces of first elements minus 1. + l2 = map(i2 + 1) + if (map(i2 + 2) <= il) then + i1 = i2 + l1 = l2 + i2 = i2 + l2 + ! if sub-sub-maps are empty, i. e. length less 4, than trara2=0 + elseif ((l1 < 4) .and. (l2 < 4)) then + trara2 = 0.0_wp + done = .true. + return + else + ! if flog2 less flog1, than ls2 first map and ls1 second map + if (map(i2 + 3) <= map(i1 + 3)) exit + call task3(done) + return + end if + end do + call task2(done) + end subroutine task1 + recursive subroutine task2(done) + logical, intent(out) :: done + done = .false. + kt = i1 i1 = i2 + i2 = kt + kt = l1 l1 = l2 - i2 = i2 + l2 - ! if sub-sub-maps are empty, i. e. length less 4, than trara2=0 - elseif ( (l1<4) .and. (l2<4) ) then - trara2 = 0.0_wp - done = .true. - return - else - ! if flog2 less flog1, than ls2 first map and ls1 second map - if ( map(i2+3)<=map(i1+3) ) exit + l2 = kt call task3(done) - return - endif - enddo - call task2(done) - end subroutine task1 - recursive subroutine task2(done) - logical,intent(out) :: done - done = .false. - kt = i1 - i1 = i2 - i2 = kt - kt = l1 - l1 = l2 - l2 = kt - call task3(done) - end subroutine task2 - recursive subroutine task3(done) - logical,intent(out) :: done - logical :: check - done = .false. - ! determine interpolate in scaled l-value - fll1 = map(i1+2) - fll2 = map(i2+2) - dfl = (fnl-fll1)/(fll2-fll1) - flog1 = map(i1+3) - flog2 = map(i2+3) - fkb1 = 0.0_wp - fkb2 = 0.0_wp - if ( l1>=4 ) then - ! b/b0 loop - check = .true. - do j2 = 4 , l2 - fincr2 = map(i2+j2) - if ( fkb2+fincr2>fnb ) then - check = .false. - exit + end subroutine task2 + recursive subroutine task3(done) + logical, intent(out) :: done + logical :: check + done = .false. + ! determine interpolate in scaled l-value + fll1 = map(i1 + 2) + fll2 = map(i2 + 2) + dfl = (fnl - fll1) / (fll2 - fll1) + flog1 = map(i1 + 3) + flog2 = map(i2 + 3) + fkb1 = 0.0_wp + fkb2 = 0.0_wp + if (l1 >= 4) then + ! b/b0 loop + check = .true. + do j2 = 4, l2 + fincr2 = map(i2 + j2) + if (fkb2 + fincr2 > fnb) then + check = .false. + exit + end if + fkb2 = fkb2 + fincr2 + flog2 = flog2 - fistep + end do + if (check) then + itime = itime + 1 + if (itime == 1) then + call task2(done) + return + end if + trara2 = 0.0_wp + done = .true. + return + end if + if (itime /= 1) then + if (j2 == 4) then + call task4(done) + return + end if + sl2 = flog2 / fkb2 + check = .true. + do j1 = 4, l1 + fincr1 = map(i1 + j1) + fkb1 = fkb1 + fincr1 + flog1 = flog1 - fistep + fkbj1 = ((flog1 / fistep) * fincr1 + fkb1) / ((fincr1 / fistep) * sl2 + 1.0_wp) + if (fkbj1 <= fkb1) then + check = .false. + exit + end if + end do + if (check) then + if (fkbj1 <= fkb2) then + trara2 = 0.0_wp + done = .true. + return + end if + end if + if (fkbj1 <= fkb2) then + fkbm = fkbj1 + (fkb2 - fkbj1) * dfl + flogm = fkbm * sl2 + flog2 = flog2 - fistep + fkb2 = fkb2 + fincr2 + sl1 = flog1 / fkb1 + sl2 = flog2 / fkb2 + call task5(done) + return + else + fkb1 = 0.0_wp + end if + end if + fkb2 = 0.0_wp end if + j2 = 4 + fincr2 = map(i2 + j2) + flog2 = map(i2 + 3) + flog1 = map(i1 + 3) + call task4(done) + end subroutine task3 + recursive subroutine task4(done) + logical, intent(out) :: done + done = .false. + flogm = flog1 + (flog2 - flog1) * dfl + fkbm = 0.0_wp fkb2 = fkb2 + fincr2 flog2 = flog2 - fistep - enddo - if (check) then - itime = itime + 1 - if ( itime==1 ) then - call task2(done) - return - endif - trara2 = 0.0_wp - done = .true. - return - end if - if ( itime/=1 ) then - if ( j2==4 ) then - call task4(done) - return - endif - sl2 = flog2/fkb2 - check = .true. - do j1 = 4 , l1 - fincr1 = map(i1+j1) - fkb1 = fkb1 + fincr1 - flog1 = flog1 - fistep - fkbj1 = ((flog1/fistep)*fincr1+fkb1)/((fincr1/fistep)*sl2+1.0_wp) - if ( fkbj1<=fkb1 ) then - check = .false. - exit - end if - enddo - if (check) then - if ( fkbj1<=fkb2 ) then - trara2 = 0.0_wp - done = .true. - return - endif - end if - if ( fkbj1<=fkb2 ) then - fkbm = fkbj1 + (fkb2-fkbj1)*dfl - flogm = fkbm*sl2 - flog2 = flog2 - fistep - fkb2 = fkb2 + fincr2 - sl1 = flog1/fkb1 - sl2 = flog2/fkb2 - call task5(done) - return + sl2 = flog2 / fkb2 + if (l1 < 4) then + fincr1 = 0.0_wp + sl1 = -900000.0_wp + call task6(done) + return else - fkb1 = 0.0_wp - endif - endif - fkb2 = 0.0_wp - endif - j2 = 4 - fincr2 = map(i2+j2) - flog2 = map(i2+3) - flog1 = map(i1+3) - call task4(done) - end subroutine task3 - recursive subroutine task4(done) - logical,intent(out) :: done - done = .false. - flogm = flog1 + (flog2-flog1)*dfl - fkbm = 0.0_wp - fkb2 = fkb2 + fincr2 - flog2 = flog2 - fistep - sl2 = flog2/fkb2 - if ( l1<4 ) then - fincr1 = 0.0_wp - sl1 = -900000.0_wp - call task6(done) - return - else - j1 = 4 - fincr1 = map(i1+j1) - fkb1 = fkb1 + fincr1 - flog1 = flog1 - fistep - sl1 = flog1/fkb1 - endif - call task5(done) - end subroutine task4 - recursive subroutine task5(done) - logical,intent(out) :: done - done = .false. - do while ( sl1>=sl2 ) - fkbj2 = ((flog2/fistep)*fincr2+fkb2)/((fincr2/fistep)*sl1+1.0_wp) - fkb = fkb1 + (fkbj2-fkb1)*dfl - flog = fkb*sl1 - if ( fkb>=fnb ) then + j1 = 4 + fincr1 = map(i1 + j1) + fkb1 = fkb1 + fincr1 + flog1 = flog1 - fistep + sl1 = flog1 / fkb1 + end if + call task5(done) + end subroutine task4 + recursive subroutine task5(done) + logical, intent(out) :: done + done = .false. + do while (sl1 >= sl2) + fkbj2 = ((flog2 / fistep) * fincr2 + fkb2) / ((fincr2 / fistep) * sl1 + 1.0_wp) + fkb = fkb1 + (fkbj2 - fkb1) * dfl + flog = fkb * sl1 + if (fkb >= fnb) then + call task7(done) + return + end if + fkbm = fkb + flogm = flog + if (j1 >= l1) then + trara2 = 0.0_wp + done = .true. + return + else + j1 = j1 + 1 + fincr1 = map(i1 + j1) + flog1 = flog1 - fistep + fkb1 = fkb1 + fincr1 + sl1 = flog1 / fkb1 + end if + end do + call task6(done) + end subroutine task5 + recursive subroutine task6(done) + logical, intent(out) :: done + done = .false. + fkbj1 = ((flog1 / fistep) * fincr1 + fkb1) / ((fincr1 / fistep) * sl2 + 1.0_wp) + fkb = fkbj1 + (fkb2 - fkbj1) * dfl + flog = fkb * sl2 + if (fkb < fnb) then + fkbm = fkb + flogm = flog + if (j2 >= l2) then + trara2 = 0.0_wp + done = .true. + return + else + j2 = j2 + 1 + fincr2 = map(i2 + j2) + flog2 = flog2 - fistep + fkb2 = fkb2 + fincr2 + sl2 = flog2 / fkb2 + call task5(done) + return + end if + end if call task7(done) - return - endif - fkbm = fkb - flogm = flog - if ( j1>=l1 ) then - trara2 = 0.0_wp - done = .true. - return - else - j1 = j1 + 1 - fincr1 = map(i1+j1) - flog1 = flog1 - fistep - fkb1 = fkb1 + fincr1 - sl1 = flog1/fkb1 - endif - enddo - call task6(done) - end subroutine task5 - recursive subroutine task6(done) - logical,intent(out) :: done - done = .false. - fkbj1 = ((flog1/fistep)*fincr1+fkb1)/((fincr1/fistep)*sl2+1.0_wp) - fkb = fkbj1 + (fkb2-fkbj1)*dfl - flog = fkb*sl2 - if ( fkb=l2 ) then - trara2 = 0.0_wp + end subroutine task6 + recursive subroutine task7(done) + logical, intent(out) :: done + if (fkb < fkbm + 1.0e-10_wp) then + trara2 = 0.0_wp + else + trara2 = flogm + (flog - flogm) * ((fnb - fkbm) / (fkb - fkbm)) + trara2 = max(trara2, 0.0_wp) + end if done = .true. - return - else - j2 = j2 + 1 - fincr2 = map(i2+j2) - flog2 = flog2 - fistep - fkb2 = fkb2 + fincr2 - sl2 = flog2/fkb2 - call task5(done) - return - endif - endif - call task7(done) - end subroutine task6 - recursive subroutine task7(done) - logical,intent(out) :: done - if ( fkb