text:
C     This File is Automatically generated by ALOHA
C     The process calculated in this file is:
C     Gamma(3,2,-1)*ProjM(-1,1) + 2*Gamma(3,2,-1)*ProjP(-1,1)
C
      SUBROUTINE FFV4_3(F1, F2, COUP, M3, W3,V3)
      IMPLICIT NONE
      COMPLEX*16 CI
      PARAMETER (CI=(0D0,1D0))
      COMPLEX*16 DENOM
      COMPLEX*16 V3(6)
      COMPLEX*16 TMP1
      REAL*8 W3
      REAL*8 P3(0:3)
      REAL*8 M3
      COMPLEX*16 F1(*)
      COMPLEX*16 F2(*)
      REAL*8 OM3
      COMPLEX*16 TMP4
      COMPLEX*16 COUP
      OM3 = 0D0
      IF (M3.NE.0D0) OM3=1D0/M3**2
      V3(1) = +F1(1)+F2(1)
      V3(2) = +F1(2)+F2(2)
      P3(0) = -DBLE(V3(1))
      P3(1) = -DBLE(V3(2))
      P3(2) = -DIMAG(V3(2))
      P3(3) = -DIMAG(V3(1))
      TMP4 = (F1(5)*(F2(3)*(P3(0)-P3(3))-F2(4)*(P3(1)+CI*(P3(2))))
     $ +F1(6)*(F2(3)*(+CI*(P3(2))-P3(1))+F2(4)*(P3(0)+P3(3))))
      TMP1 = (F1(3)*(F2(5)*(P3(0)+P3(3))+F2(6)*(P3(1)+CI*(P3(2))))
     $ +F1(4)*(F2(5)*(P3(1)-CI*(P3(2)))+F2(6)*(P3(0)-P3(3))))
      DENOM = COUP/(P3(0)**2-P3(1)**2-P3(2)**2-P3(3)**2 - M3 * (M3 -CI
     $ * W3))
      V3(3)= DENOM*(-2D0 * CI)*(OM3*-1D0/2D0 * P3(0)*(TMP1+2D0*(TMP4))
     $ +(+1D0/2D0*(F1(3)*F2(5)+F1(4)*F2(6))+F1(5)*F2(3)+F1(6)*F2(4)))
      V3(4)= DENOM*(-2D0 * CI)*(OM3*-1D0/2D0 * P3(1)*(TMP1+2D0*(TMP4))
     $ +(-1D0/2D0*(F1(3)*F2(6)+F1(4)*F2(5))+F1(5)*F2(4)+F1(6)*F2(3)))
      V3(5)= DENOM*2D0 * CI*(OM3*1D0/2D0 * P3(2)*(TMP1+2D0*(TMP4))+(
     $ +1D0/2D0 * CI*(F1(3)*F2(6))-1D0/2D0 * CI*(F1(4)*F2(5))-CI*(F1(5)
     $ *F2(4))+CI*(F1(6)*F2(3))))
      V3(6)= DENOM*2D0 * CI*(OM3*1D0/2D0 * P3(3)*(TMP1+2D0*(TMP4))+(
     $ +1D0/2D0*(F1(3)*F2(5))-1D0/2D0*(F1(4)*F2(6))-F1(5)*F2(3)+F1(6)
     $ *F2(4)))
      END

max_stars_repo_name: madnklo/madnklo

max_issues_repo_name: madnklo/madnklo

max_forks_repo_name: madnklo/madnklo


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!

subroutine psi_renum_index(iperm,idx,info)
  use psi_mod, psi_protect_name =>  psi_renum_index
  use psb_serial_mod
  implicit none

  integer(psb_ipk_), intent(out)   :: info
  integer(psb_ipk_), intent(in)    :: iperm(:)
  integer(psb_ipk_), intent(inout) :: idx(:)

  integer(psb_ipk_) :: i,j,k,nh

  i=1
  k=idx(i)
  do while (k /= -1)
    i = i+1
    nh = idx(i)
    do j = i+1, i+nh
      idx(j) = iperm(idx(j))
    enddo
    i  = i + nh + 1
    nh = idx(i)
    do j = i+1, i+nh
      idx(j) = iperm(idx(j))
    enddo
    i = i + nh + 1
    k = idx(i)
  enddo

end subroutine psi_renum_index

subroutine psi_i_cnv_dsc(halo_in,ovrlap_in,ext_in,cdesc, info, mold)

  use psi_mod, psi_protect_name =>  psi_i_cnv_dsc
  use psb_realloc_mod
  implicit none

  !     ....scalars parameters....
  integer(psb_ipk_), intent(in)      :: halo_in(:), ovrlap_in(:),ext_in(:)
  type(psb_desc_type), intent(inout) :: cdesc
  integer(psb_ipk_), intent(out)     :: info
  class(psb_i_base_vect_type), optional, intent(in) :: mold

  !     ....local scalars....
  integer(psb_ipk_) :: np,me
  integer(psb_ipk_) :: ictxt, err_act,nxch,nsnd,nrcv,j,k
  !     ...local array...
  integer(psb_ipk_), allocatable  :: idx_out(:), tmp_mst_idx(:)

  !     ...parameters
  integer(psb_ipk_) :: debug_level, debug_unit
  logical, parameter :: debug=.false.
  character(len=20)  :: name

  name='psi_cnv_desc'
  call psb_get_erraction(err_act)
  debug_level = psb_get_debug_level()
  debug_unit  = psb_get_debug_unit()

  info = psb_success_
  ictxt = cdesc%get_context()

  call psb_info(ictxt,me,np)
  if (np == -1) then
    info = psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  ! first the halo index
  if (debug_level>0) write(debug_unit,*) me,'Calling crea_index on halo',&
       & size(halo_in)
  call psi_crea_index(cdesc,halo_in, idx_out,nxch,nsnd,nrcv,info)
  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_crea_index')
    goto 9999
  end if
  call psb_move_alloc(idx_out,cdesc%halo_index,info)

  if (debug_level>0) write(debug_unit,*) me,'Done crea_index on halo'
  if (debug_level>0) write(debug_unit,*) me,'Calling crea_index on ext'

  ! then ext index
  if (debug_level>0) write(debug_unit,*) me,'Calling crea_index on ext'
  call psi_crea_index(cdesc,ext_in, idx_out,nxch,nsnd,nrcv,info)
  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_crea_index')
    goto 9999
  end if
  call psb_move_alloc(idx_out,cdesc%ext_index,info)

  if (debug_level>0) write(debug_unit,*) me,'Done crea_index on ext'
  if (debug_level>0) write(debug_unit,*) me,'Calling crea_index on ovrlap'

  ! then the overlap index
  call psi_crea_index(cdesc,ovrlap_in, idx_out,nxch,nsnd,nrcv,info)
  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_crea_index')
    goto 9999
  end if
  call psb_move_alloc(idx_out,cdesc%ovrlap_index,info)
  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_move_alloc')
    goto 9999
  end if

  ! next  ovrlap_elem
  if (debug_level>0) write(debug_unit,*) me,'Calling crea_ovr_elem'
  call psi_crea_ovr_elem(me,cdesc%ovrlap_index,cdesc%ovrlap_elem,info)
  if (debug_level>0) write(debug_unit,*) me,'Done crea_ovr_elem'
  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_crea_ovr_elem')
    goto 9999
  end if
  ! Extract ovr_mst_idx from ovrlap_elem
  if (debug_level>0) write(debug_unit,*) me,'Calling bld_ovr_mst'
  call psi_bld_ovr_mst(me,cdesc%ovrlap_elem,tmp_mst_idx,info)
  if (info == psb_success_) call psi_crea_index(cdesc,&
       & tmp_mst_idx,idx_out,nxch,nsnd,nrcv,info)
  if (debug_level>0) write(debug_unit,*) me,'Done crea_indx'
  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_bld_ovr_mst')
    goto 9999
  end if
  call psb_move_alloc(idx_out,cdesc%ovr_mst_idx,info)
  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_move_alloc')
    goto 9999
  end if

  ! finally bnd_elem
  call psi_crea_bnd_elem(idx_out,cdesc,info)
  if (info == psb_success_) call psb_move_alloc(idx_out,cdesc%bnd_elem,info)

  call cdesc%v_halo_index%bld(cdesc%halo_index,mold=mold)
  call cdesc%v_ext_index%bld(cdesc%ext_index,mold=mold)
  call cdesc%v_ovrlap_index%bld(cdesc%ovrlap_index,mold=mold)
  call cdesc%v_ovr_mst_idx%bld(cdesc%ovr_mst_idx,mold=mold)

  if (info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_crea_bnd_elem')
    goto 9999
  end if
  if (debug_level>0) write(debug_unit,*) me,'Done crea_bnd_elem'

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(ictxt,err_act)

  return

end subroutine psi_i_cnv_dsc

subroutine psi_i_inner_cnvs(x,hashmask,hashv,glb_lc)
  use psi_mod, psi_protect_name => psi_i_inner_cnvs

  integer(psb_ipk_), intent(in)    :: hashmask,hashv(0:),glb_lc(:,:)
  integer(psb_ipk_), intent(inout) :: x

  integer(psb_ipk_) :: i, ih, key, idx,nh,tmp,lb,ub,lm
  !
  ! When a large descriptor is assembled the indices
  ! are kept in a (hashed) list of ordered lists.
  ! Thus we first hash the index, then we do a binary search on the
  ! ordered sublist. The hashing is based on the low-order bits
  ! for a width of psb_hash_bits
  !

  key = x
  ih  = iand(key,hashmask)
  idx = hashv(ih)
  nh  = hashv(ih+1) - hashv(ih)
  if (nh > 0) then
    tmp = -1
    lb = idx
    ub = idx+nh-1
    do
      if (lb>ub) exit
      lm = (lb+ub)/2
      if (key == glb_lc(lm,1)) then
        tmp = lm
        exit
      else if (key<glb_lc(lm,1)) then
        ub = lm - 1
      else
        lb = lm + 1
      end if
    end do
  else
    tmp = -1
  end if
  if (tmp > 0) then
    x = glb_lc(tmp,2)
  else
    x = tmp
  end if
end subroutine psi_i_inner_cnvs

subroutine psi_i_inner_cnvs2(x,y,hashmask,hashv,glb_lc)
  use psi_mod, psi_protect_name =>  psi_i_inner_cnvs2
  integer(psb_ipk_), intent(in)  :: hashmask,hashv(0:),glb_lc(:,:)
  integer(psb_ipk_), intent(in)  :: x
  integer(psb_ipk_), intent(out) :: y

  integer(psb_ipk_) :: i, ih, key, idx,nh,tmp,lb,ub,lm
  !
  ! When a large descriptor is assembled the indices
  ! are kept in a (hashed) list of ordered lists.
  ! Thus we first hash the index, then we do a binary search on the
  ! ordered sublist. The hashing is based on the low-order bits
  ! for a width of psb_hash_bits
  !

  key = x
  ih  = iand(key,hashmask)
  idx = hashv(ih)
  nh  = hashv(ih+1) - hashv(ih)
  if (nh > 0) then
    tmp = -1
    lb = idx
    ub = idx+nh-1
    do
      if (lb>ub) exit
      lm = (lb+ub)/2
      if (key == glb_lc(lm,1)) then
        tmp = lm
        exit
      else if (key<glb_lc(lm,1)) then
        ub = lm - 1
      else
        lb = lm + 1
      end if
    end do
  else
    tmp = -1
  end if
  if (tmp > 0) then
    y = glb_lc(tmp,2)
  else
    y = tmp
  end if
end subroutine psi_i_inner_cnvs2

subroutine psi_i_inner_cnv1(n,x,hashmask,hashv,glb_lc,mask)
  use psi_mod, psi_protect_name =>  psi_i_inner_cnv1
  integer(psb_ipk_), intent(in)    :: n,hashmask,hashv(0:),glb_lc(:,:)
  logical, intent(in), optional    :: mask(:)
  integer(psb_ipk_), intent(inout) :: x(:)

  integer(psb_ipk_) :: i, ih, key, idx,nh,tmp,lb,ub,lm
  !
  ! When a large descriptor is assembled the indices
  ! are kept in a (hashed) list of ordered lists.
  ! Thus we first hash the index, then we do a binary search on the
  ! ordered sublist. The hashing is based on the low-order bits
  ! for a width of psb_hash_bits
  !
  if (present(mask)) then
    do i=1, n
      if (mask(i)) then
        key = x(i)
        ih  = iand(key,hashmask)
        idx = hashv(ih)
        nh  = hashv(ih+1) - hashv(ih)
        if (nh > 0) then
          tmp = -1
          lb = idx
          ub = idx+nh-1
          do
            if (lb>ub) exit
            lm = (lb+ub)/2
            if (key == glb_lc(lm,1)) then
              tmp = lm
              exit
            else if (key<glb_lc(lm,1)) then
              ub = lm - 1
            else
              lb = lm + 1
            end if
          end do
        else
          tmp = -1
        end if
        if (tmp > 0) then
          x(i) = glb_lc(tmp,2)
        else
          x(i) = tmp
        end if
      end if
    end do
  else
    do i=1, n
      key = x(i)
      ih  = iand(key,hashmask)
      idx = hashv(ih)
      nh  = hashv(ih+1) - hashv(ih)
      if (nh > 0) then
        tmp = -1
        lb = idx
        ub = idx+nh-1
        do
          if (lb>ub) exit
          lm = (lb+ub)/2
          if (key == glb_lc(lm,1)) then
            tmp = lm
            exit
          else if (key<glb_lc(lm,1)) then
            ub = lm - 1
          else
            lb = lm + 1
          end if
        end do
      else
        tmp = -1
      end if
      if (tmp > 0) then
        x(i) = glb_lc(tmp,2)
      else
        x(i) = tmp
      end if
    end do
  end if
end subroutine psi_i_inner_cnv1

subroutine psi_i_inner_cnv2(n,x,y,hashmask,hashv,glb_lc,mask)
  use psi_mod, psi_protect_name =>  psi_i_inner_cnv2
  integer(psb_ipk_), intent(in)  :: n, hashmask,hashv(0:),glb_lc(:,:)
  logical, intent(in),optional  :: mask(:)
  integer(psb_ipk_), intent(in)  :: x(:)
  integer(psb_ipk_), intent(out) :: y(:)

  integer(psb_ipk_) :: i, ih, key, idx,nh,tmp,lb,ub,lm
  !
  ! When a large descriptor is assembled the indices
  ! are kept in a (hashed) list of ordered lists.
  ! Thus we first hash the index, then we do a binary search on the
  ! ordered sublist. The hashing is based on the low-order bits
  ! for a width of psb_hash_bits
  !
  if (present(mask)) then
    do i=1, n
      if (mask(i)) then
        key = x(i)
        ih  = iand(key,hashmask)
        if (ih > ubound(hashv,1) ) then
          write(psb_err_unit,*) ' In inner cnv: ',ih,ubound(hashv)
        end if
        idx = hashv(ih)
        nh  = hashv(ih+1) - hashv(ih)
        if (nh > 0) then
          tmp = -1
          lb = idx
          ub = idx+nh-1
          do
            if (lb>ub) exit
            lm = (lb+ub)/2
            if (key == glb_lc(lm,1)) then
              tmp = lm
              exit
            else if (key<glb_lc(lm,1)) then
              ub = lm - 1
            else
              lb = lm + 1
            end if
          end do
        else
          tmp = -1
        end if
        if (tmp > 0) then
          y(i) = glb_lc(tmp,2)
        else
          y(i) = tmp
        end if
      end if
    end do
  else
    do i=1, n
      key = x(i)
      ih  = iand(key,hashmask)
      if (ih > ubound(hashv,1) ) then
        write(psb_err_unit,*) ' In inner cnv: ',ih,ubound(hashv)
      end if
      idx = hashv(ih)
      nh  = hashv(ih+1) - hashv(ih)
      if (nh > 0) then
        tmp = -1
        lb = idx
        ub = idx+nh-1
        do
          if (lb>ub) exit
          lm = (lb+ub)/2
          if (key == glb_lc(lm,1)) then
            tmp = lm
            exit
          else if (key<glb_lc(lm,1)) then
            ub = lm - 1
          else
            lb = lm + 1
          end if
        end do
      else
        tmp = -1
      end if
      if (tmp > 0) then
        y(i) = glb_lc(tmp,2)
      else
        y(i) = tmp
      end if
    end do
  end if
end subroutine psi_i_inner_cnv2

subroutine psi_i_bld_ovr_mst(me,ovrlap_elem,mst_idx,info)
  use psi_mod, psi_protect_name =>  psi_i_bld_ovr_mst

  use psb_realloc_mod
  implicit none

  !     ....scalars parameters....
  integer(psb_ipk_), intent(in)               :: me, ovrlap_elem(:,:)
  integer(psb_ipk_), allocatable, intent(out) :: mst_idx(:)
  integer(psb_ipk_), intent(out)              :: info

  integer(psb_ipk_) :: i, j, proc, nov,isz, ip, err_act, idx
  character(len=20)  :: name

  name='psi_bld_ovr_mst'
  call psb_get_erraction(err_act)

  nov = size(ovrlap_elem,1)
  isz = 3*nov+1
  call psb_realloc(isz,mst_idx,info)
  if (info /= psb_success_) then
    call psb_errpush(psb_err_internal_error_,name,a_err='reallocate')
    goto 9999
  end if
  mst_idx = -1
  j = 1
  do i=1, nov
    proc = ovrlap_elem(i,3)
    if (me /= proc) then
      idx = ovrlap_elem(i,1)
      mst_idx(j+0) = proc
      mst_idx(j+1) = 1
      mst_idx(j+2) = idx
      j = j + 3
    end if
  end do
  mst_idx(j) = -1

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return

end subroutine psi_i_bld_ovr_mst


text:
MODULE libtest

INTERFACE
        SUBROUTINE test_exit(exit_status) BIND(C, NAME='exit')
                USE, INTRINSIC :: ISO_C_BINDING
                IMPLICIT NONE
                INTEGER (C_INT), INTENT(IN), VALUE :: exit_status
        ENDSUBROUTINE test_exit
ENDINTERFACE

CONTAINS

SUBROUTINE ASSERT(condition)
        USE, INTRINSIC :: ISO_FORTRAN_ENV
        IMPLICIT NONE
        LOGICAL, INTENT(IN) :: condition
        INTEGER, SAVE :: counter = 1
        WRITE(output_unit, '(A, I5)') 'assert call ', counter
        counter = counter + 1
        IF (.NOT. condition) CALL test_exit(1)
ENDSUBROUTINE ASSERT
ENDMODULE libtest

max_stars_repo_name: mschollmeier/xraylib

max_issues_repo_name: mschollmeier/xraylib

max_forks_repo_name: mschollmeier/xraylib


text:
      double precision function dqwgtf(x,omega,p2,p3,p4,integr)
c***begin prologue  dqwgtf
c***refer to   dqk15w
c***routines called  (none)
c***revision date 810101   (yymmdd)
c***keywords  cos or sin in weight function
c***author  piessens,robert, appl. math. & progr. div. - k.u.leuven
c           de doncker,elise,appl. math. * progr. div. - k.u.leuven
c***end prologue  dqwgtf
c
      double precision dcos,dsin,omega,omx,p2,p3,p4,x
      integer integr
c***first executable statement  dqwgtf
      omx = omega*x
      go to(10,20),integr
   10 dqwgtf = dcos(omx)
      go to 30
   20 dqwgtf = dsin(omx)
   30 return
      end

ext: f

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy


text:
      SUBROUTINE E03118 (USEPRM)

C  E03118 tests the handling of error 118.

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /ERRINF/ ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR,
     1                USRERR, ERRSAV,      FUNSAV,      FILSAV,
     2                EFCNT, EFID
      INTEGER         ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR(10),
     1                USRERR, ERRSAV(200), FUNSAV(200), FILSAV(200),
     2                EFCNT, EFID(100)
      COMMON /ERRCHR/ CURCON,     ERRSRS,    ERRMRK,    ERFLNM,
     1                CONTAB
      CHARACTER       CURCON*200, ERRSRS*40, ERRMRK*20, ERFLNM*80,
     1                CONTAB(40)*150

C  type of returned values
      INTEGER        PSET,      PREALI
      PARAMETER     (PSET=0,    PREALI=1)

C colour model
      INTEGER    PRGB,     PCIE,     PHSV,     PHLS
      PARAMETER (PRGB = 1, PCIE = 2, PHSV = 3, PHLS = 4)

      INTEGER   USEPRM, NCLMOD, IDUM1, IDUM2, LCLMOD, IX
      INTEGER   NUMCP2, ISIZ, SPECWT
      REAL      SPEC1(10), SPEC2(10), SPEC(10)
      LOGICAL   STREQ,  RAREQ

      CURCON = 'one of the components of the colour ' //
     1         'specification is out of range'
      CALL SETVS ('118', EXPERR, EXPSIZ)
      ERRSRS = '17'
      CALL ESETUP (USEPRM)

      CALL POPWK (WKID, CONID, WTYPE)
      CALL PQWKC (WKID, ERRIND, CONID, SPECWT)
      CALL CHKINQ ('pqwkc', ERRIND)

      CALL PQCMDF (SPECWT,0,ERRIND,NCLMOD,IDUM1,IDUM2)
      CALL CHKINQ ('pqcmdf', ERRIND)

      DO 150 IX = 1, NCLMOD

         CALL PQCMDF (SPECWT, IX, ERRIND, IDUM1, LCLMOD,IDUM2)
         CALL CHKINQ ('pqcmdf', ERRIND)
         IF (LCLMOD.LT.PRGB .OR. LCLMOD.GT.PHLS) GOTO 150

         CALL PSCMD (WKID, LCLMOD)
         CALL SETRVS ('0.5,0.5,0.5', SPEC1, ISIZ)
         CALL PSCR (WKID, 2, 3, SPEC1)

         CALL SETRVS ('-0.5, 0., 1.', SPEC, ISIZ)
         CALL RSCR (WKID, 2, 3, SPEC)

         CALL SETRVS ('0., 1.01, 0.5', SPEC, ISIZ)
         CALL RSCR (WKID, 2, 3, SPEC)

         CALL SETRVS ('0., 0.5, -0.01', SPEC, ISIZ)
         CALL RSCR (WKID, 2, 3, SPEC)

         CALL PQCR (WKID,2,10,PSET,ERRIND,NUMCP2,SPEC2)
         CALL CHKINQ ('pqcr', ERRIND)
         CALL TSTIGN (STREQ('OO**')  .AND.
     1                NUMCP2.EQ.3    .AND.
     2                RAREQ(3,SPEC2,SPEC1,0.01,0.01))

150   CONTINUE

      CALL PCLWK (WKID)

      CALL ENDERR

      END


text:
	module sub
	contains

c	This subroutine is used to calculate p and grad_p which use the data of
c	latitude and longtitude
c	x means lon,y means lat,N_max means Nmax here
	subroutine point(N_max,x,y,p,grad_p)
	implicit none
	real*8 x
	real*8 y
	real*8 p(N_max*(N_max+3)/2+1)
	real*8 grad_p(N_max*(N_max+3)/2+1)
	integer N_max,n,m,i

	p(1)=1
	p(2)=cos(y)
	p(3)=sin(y)
	grad_p(1)=0
	grad_p(2)=(cos(y)*p(2)-p(1))/sin(y)
	grad_p(3)=cos(y)*p(3)/sin(y)
	if (N_max>=2) then
	i=4
	do n=2,N_max
	do m=0,n
	if(m==n)then
	p(i)=p(i-n-1)*sin(y)*((2*m-1)/(2*m))**0.5
	grad_p(i)=n*cos(y)*p(i-n-1)*((2*m-1)/(2*m))**0.5
	else
	p(i)=((2*n-1)*cos(y)*p(i-n)-((n-1)**2-m**2)**
     &	0.5*p(i-n-n+1))/(n**2-m**2)**0.5
	grad_p(i)=(n*cos(y)*p(i)-(n**2-m**2)**
     &	0.5*p(i-n))/sin(y)
	end if
	i=i+1
	end do
	end do
	end if
	end subroutine point

c	得到观测方程的系数矩阵LA，j表示计算LA第几行，即第j个观测点的系数行。
	subroutine getLAofx(lat,lon,r,N_max,LA,data_N,j)
	implicit none
	integer N_max,j,data_N
	real*8 lat,lon,r,LA(data_N,N_max*(N_max+2))
	real*8 p(N_max*(N_max+3)/2+1),grad_p(N_max*(N_max+3)/2+1)
	integer i,n,m,k
	call point(N_max,lon,lat,p,grad_p)
c	ux=ux+(g(i)*cos(m*x)+h(i)*sin(m*x))*grad_p(i+1)
	i=1
	k=1
	do n=1,N_max
	do m=0,n
	LA(j,i)=cos(m*lon)*grad_p(i+1)
	if(m>0) then
	LA(j,k+N_max*(N_max+3)/2)=sin(m*lon)*grad_p(i+1)
	k=k+1
	end if
	i=i+1
	end do
	end do
	call igrf11syn(
	end subroutine

      subroutine igrf11syn (isv,date,itype,alt,colat,elong,x,y,z,f)
c
c     This is a synthesis routine for the 11th generation IGRF as agreed
c     in December 2009 by IAGA Working Group V-MOD. It is valid 1900.0 to
c     2015.0 inclusive. Values for dates from 1945.0 to 2005.0 inclusive are
c     definitive, otherwise they are non-definitive.
c   INPUT
c     isv   = 0 if main-field values are required
c     isv   = 1 if secular variation values are required
c     date  = year A.D. Must be greater than or equal to 1900.0 and
c             less than or equal to 2020.0. Warning message is given
c             for dates greater than 2015.0. Must be double precision.
c     itype = 1 if geodetic (spheroid)
c     itype = 2 if geocentric (sphere)
c     alt   = height in km above sea level if itype = 1
c           = distance from centre of Earth in km if itype = 2 (>3485 km)
c     colat = colatitude (0-180)
c     elong = east-longitude (0-360)
c     alt, colat and elong must be double precision.
c   OUTPUT
c     x     = north component (nT) if isv = 0, nT/year if isv = 1
c     y     = east component (nT) if isv = 0, nT/year if isv = 1
c     z     = vertical component (nT) if isv = 0, nT/year if isv = 1
c     f     = total intensity (nT) if isv = 0, rubbish if isv = 1
c
c     To get the other geomagnetic elements (D, I, H and secular
c     variations dD, dH, dI and dF) use routines ptoc and ptocsv.
c
c     Adapted from 8th generation version to include new maximum degree for
c     main-field models for 2000.0 and onwards and use WGS84 spheroid instead
c     of International Astronomical Union 1966 spheroid as recommended by IAGA
c     in July 2003. Reference radius remains as 6371.2 km - it is NOT the mean
c     radius (= 6371.0 km) but 6371.2 km is what is used in determining the
c     coefficients. Adaptation by Susan Macmillan, August 2003 (for
c     9th generation), December 2004 & December 2009.
c
c     Coefficients at 1995.0 incorrectly rounded (rounded up instead of
c     to even) included as these are the coefficients published in Excel
c     spreadsheet July 2005.
c
      implicit double precision (a-h,o-z)
      dimension gh(3256),g0(120),g1(120),g2(120),g3(120),g4(120),
     1          g5(120),g6(120),g7(120),g8(120),g9(120),ga(120),
     2          gb(120),gc(120),gd(120),ge(120),gf(120),gg(120),
     3          gi(120),gj(120),gk(195),gl(195),gm(195),gp(195),
     4          gq(195),
     5          p(105),q(105),cl(13),sl(13)
      equivalence (g0,gh(1)),(g1,gh(121)),(g2,gh(241)),(g3,gh(361)),
     1            (g4,gh(481)),(g5,gh(601)),(g6,gh(721)),(g7,gh(841)),
     2            (g8,gh(961)),(g9,gh(1081)),(ga,gh(1201)),
     3            (gb,gh(1321)),(gc,gh(1441)),(gd,gh(1561)),
     4            (ge,gh(1681)),(gf,gh(1801)),(gg,gh(1921)),
     5            (gi,gh(2041)),(gj,gh(2161)),(gk,gh(2281)),
     6            (gl,gh(2476)),(gm,gh(2671)),(gp,gh(2866)),
     7            (gq,gh(3061))
c
      data g0/ -31543.,-2298., 5922., -677., 2905.,-1061.,  924., 1121., 1900
     1           1022.,-1469., -330., 1256.,    3.,  572.,  523.,  876., 1900
     2            628.,  195.,  660.,  -69., -361., -210.,  134.,  -75., 1900
     3           -184.,  328., -210.,  264.,   53.,    5.,  -33.,  -86., 1900
     4           -124.,  -16.,    3.,   63.,   61.,   -9.,  -11.,   83., 1900
     5           -217.,    2.,  -58.,  -35.,   59.,   36.,  -90.,  -69., 1900
     6             70.,  -55.,  -45.,    0.,  -13.,   34.,  -10.,  -41., 1900
     7             -1.,  -21.,   28.,   18.,  -12.,    6.,  -22.,   11., 1900
     8              8.,    8.,   -4.,  -14.,   -9.,    7.,    1.,  -13., 1900
     9              2.,    5.,   -9.,   16.,    5.,   -5.,    8.,  -18., 1900
     a              8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., 1900
     b             -3.,    1.,   -2.,   -2.,    8.,    2.,   10.,   -1., 1900
     c             -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., 1900
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1900
     e              0.,   -2.,    2.,    4.,    2.,    0.,    0.,   -6./ 1900
      data g1/ -31464.,-2298., 5909., -728., 2928.,-1086., 1041., 1065., 1905
     1           1037.,-1494., -357., 1239.,   34.,  635.,  480.,  880., 1905
     2            643.,  203.,  653.,  -77., -380., -201.,  146.,  -65., 1905
     3           -192.,  328., -193.,  259.,   56.,   -1.,  -32.,  -93., 1905
     4           -125.,  -26.,   11.,   62.,   60.,   -7.,  -11.,   86., 1905
     5           -221.,    4.,  -57.,  -32.,   57.,   32.,  -92.,  -67., 1905
     6             70.,  -54.,  -46.,    0.,  -14.,   33.,  -11.,  -41., 1905
     7              0.,  -20.,   28.,   18.,  -12.,    6.,  -22.,   11., 1905
     8              8.,    8.,   -4.,  -15.,   -9.,    7.,    1.,  -13., 1905
     9              2.,    5.,   -8.,   16.,    5.,   -5.,    8.,  -18., 1905
     a              8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., 1905
     b             -3.,    1.,   -2.,   -2.,    8.,    2.,   10.,    0., 1905
     c             -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., 1905
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1905
     e              0.,   -2.,    2.,    4.,    2.,    0.,    0.,   -6./ 1905
      data g2/ -31354.,-2297., 5898., -769., 2948.,-1128., 1176., 1000., 1910
     1           1058.,-1524., -389., 1223.,   62.,  705.,  425.,  884., 1910
     2            660.,  211.,  644.,  -90., -400., -189.,  160.,  -55., 1910
     3           -201.,  327., -172.,  253.,   57.,   -9.,  -33., -102., 1910
     4           -126.,  -38.,   21.,   62.,   58.,   -5.,  -11.,   89., 1910
     5           -224.,    5.,  -54.,  -29.,   54.,   28.,  -95.,  -65., 1910
     6             71.,  -54.,  -47.,    1.,  -14.,   32.,  -12.,  -40., 1910
     7              1.,  -19.,   28.,   18.,  -13.,    6.,  -22.,   11., 1910
     8              8.,    8.,   -4.,  -15.,   -9.,    6.,    1.,  -13., 1910
     9              2.,    5.,   -8.,   16.,    5.,   -5.,    8.,  -18., 1910
     a              8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., 1910
     b             -3.,    1.,   -2.,   -2.,    8.,    2.,   10.,    0., 1910
     c             -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., 1910
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1910
     e              0.,   -2.,    2.,    4.,    2.,    0.,    0.,   -6./ 1910
      data g3/ -31212.,-2306., 5875., -802., 2956.,-1191., 1309.,  917., 1915
     1           1084.,-1559., -421., 1212.,   84.,  778.,  360.,  887., 1915
     2            678.,  218.,  631., -109., -416., -173.,  178.,  -51., 1915
     3           -211.,  327., -148.,  245.,   58.,  -16.,  -34., -111., 1915
     4           -126.,  -51.,   32.,   61.,   57.,   -2.,  -10.,   93., 1915
     5           -228.,    8.,  -51.,  -26.,   49.,   23.,  -98.,  -62., 1915
     6             72.,  -54.,  -48.,    2.,  -14.,   31.,  -12.,  -38., 1915
     7              2.,  -18.,   28.,   19.,  -15.,    6.,  -22.,   11., 1915
     8              8.,    8.,   -4.,  -15.,   -9.,    6.,    2.,  -13., 1915
     9              3.,    5.,   -8.,   16.,    6.,   -5.,    8.,  -18., 1915
     a              8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., 1915
     b             -3.,    1.,   -2.,   -2.,    8.,    2.,   10.,    0., 1915
     c             -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., 1915
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1915
     e              0.,   -2.,    1.,    4.,    2.,    0.,    0.,   -6./ 1915
      data g4/ -31060.,-2317., 5845., -839., 2959.,-1259., 1407.,  823., 1920
     1           1111.,-1600., -445., 1205.,  103.,  839.,  293.,  889., 1920
     2            695.,  220.,  616., -134., -424., -153.,  199.,  -57., 1920
     3           -221.,  326., -122.,  236.,   58.,  -23.,  -38., -119., 1920
     4           -125.,  -62.,   43.,   61.,   55.,    0.,  -10.,   96., 1920
     5           -233.,   11.,  -46.,  -22.,   44.,   18., -101.,  -57., 1920
     6             73.,  -54.,  -49.,    2.,  -14.,   29.,  -13.,  -37., 1920
     7              4.,  -16.,   28.,   19.,  -16.,    6.,  -22.,   11., 1920
     8              7.,    8.,   -3.,  -15.,   -9.,    6.,    2.,  -14., 1920
     9              4.,    5.,   -7.,   17.,    6.,   -5.,    8.,  -19., 1920
     a              8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., 1920
     b             -3.,    1.,   -2.,   -2.,    9.,    2.,   10.,    0., 1920
     c             -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., 1920
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1920
     e              0.,   -2.,    1.,    4.,    3.,    0.,    0.,   -6./ 1920
      data g5/ -30926.,-2318., 5817., -893., 2969.,-1334., 1471.,  728., 1925
     1           1140.,-1645., -462., 1202.,  119.,  881.,  229.,  891., 1925
     2            711.,  216.,  601., -163., -426., -130.,  217.,  -70., 1925
     3           -230.,  326.,  -96.,  226.,   58.,  -28.,  -44., -125., 1925
     4           -122.,  -69.,   51.,   61.,   54.,    3.,   -9.,   99., 1925
     5           -238.,   14.,  -40.,  -18.,   39.,   13., -103.,  -52., 1925
     6             73.,  -54.,  -50.,    3.,  -14.,   27.,  -14.,  -35., 1925
     7              5.,  -14.,   29.,   19.,  -17.,    6.,  -21.,   11., 1925
     8              7.,    8.,   -3.,  -15.,   -9.,    6.,    2.,  -14., 1925
     9              4.,    5.,   -7.,   17.,    7.,   -5.,    8.,  -19., 1925
     a              8.,   10.,  -20.,    1.,   14.,  -11.,    5.,   12., 1925
     b             -3.,    1.,   -2.,   -2.,    9.,    2.,   10.,    0., 1925
     c             -2.,   -1.,    2.,   -3.,   -4.,    2.,    2.,    1., 1925
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1925
     e              0.,   -2.,    1.,    4.,    3.,    0.,    0.,   -6./ 1925
      data g6/ -30805.,-2316., 5808., -951., 2980.,-1424., 1517.,  644., 1930
     1           1172.,-1692., -480., 1205.,  133.,  907.,  166.,  896., 1930
     2            727.,  205.,  584., -195., -422., -109.,  234.,  -90., 1930
     3           -237.,  327.,  -72.,  218.,   60.,  -32.,  -53., -131., 1930
     4           -118.,  -74.,   58.,   60.,   53.,    4.,   -9.,  102., 1930
     5           -242.,   19.,  -32.,  -16.,   32.,    8., -104.,  -46., 1930
     6             74.,  -54.,  -51.,    4.,  -15.,   25.,  -14.,  -34., 1930
     7              6.,  -12.,   29.,   18.,  -18.,    6.,  -20.,   11., 1930
     8              7.,    8.,   -3.,  -15.,   -9.,    5.,    2.,  -14., 1930
     9              5.,    5.,   -6.,   18.,    8.,   -5.,    8.,  -19., 1930
     a              8.,   10.,  -20.,    1.,   14.,  -12.,    5.,   12., 1930
     b             -3.,    1.,   -2.,   -2.,    9.,    3.,   10.,    0., 1930
     c             -2.,   -2.,    2.,   -3.,   -4.,    2.,    2.,    1., 1930
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1930
     e              0.,   -2.,    1.,    4.,    3.,    0.,    0.,   -6./ 1930
      data g7/ -30715.,-2306., 5812.,-1018., 2984.,-1520., 1550.,  586., 1935
     1           1206.,-1740., -494., 1215.,  146.,  918.,  101.,  903., 1935
     2            744.,  188.,  565., -226., -415.,  -90.,  249., -114., 1935
     3           -241.,  329.,  -51.,  211.,   64.,  -33.,  -64., -136., 1935
     4           -115.,  -76.,   64.,   59.,   53.,    4.,   -8.,  104., 1935
     5           -246.,   25.,  -25.,  -15.,   25.,    4., -106.,  -40., 1935
     6             74.,  -53.,  -52.,    4.,  -17.,   23.,  -14.,  -33., 1935
     7              7.,  -11.,   29.,   18.,  -19.,    6.,  -19.,   11., 1935
     8              7.,    8.,   -3.,  -15.,   -9.,    5.,    1.,  -15., 1935
     9              6.,    5.,   -6.,   18.,    8.,   -5.,    7.,  -19., 1935
     a              8.,   10.,  -20.,    1.,   15.,  -12.,    5.,   11., 1935
     b             -3.,    1.,   -3.,   -2.,    9.,    3.,   11.,    0., 1935
     c             -2.,   -2.,    2.,   -3.,   -4.,    2.,    2.,    1., 1935
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1935
     e              0.,   -1.,    2.,    4.,    3.,    0.,    0.,   -6./ 1935
      data g8/ -30654.,-2292., 5821.,-1106., 2981.,-1614., 1566.,  528., 1940
     1           1240.,-1790., -499., 1232.,  163.,  916.,   43.,  914., 1940
     2            762.,  169.,  550., -252., -405.,  -72.,  265., -141., 1940
     3           -241.,  334.,  -33.,  208.,   71.,  -33.,  -75., -141., 1940
     4           -113.,  -76.,   69.,   57.,   54.,    4.,   -7.,  105., 1940
     5           -249.,   33.,  -18.,  -15.,   18.,    0., -107.,  -33., 1940
     6             74.,  -53.,  -52.,    4.,  -18.,   20.,  -14.,  -31., 1940
     7              7.,   -9.,   29.,   17.,  -20.,    5.,  -19.,   11., 1940
     8              7.,    8.,   -3.,  -14.,  -10.,    5.,    1.,  -15., 1940
     9              6.,    5.,   -5.,   19.,    9.,   -5.,    7.,  -19., 1940
     a              8.,   10.,  -21.,    1.,   15.,  -12.,    5.,   11., 1940
     b             -3.,    1.,   -3.,   -2.,    9.,    3.,   11.,    1., 1940
     c             -2.,   -2.,    2.,   -3.,   -4.,    2.,    2.,    1., 1940
     d             -5.,    2.,   -2.,    6.,    6.,   -4.,    4.,    0., 1940
     e              0.,   -1.,    2.,    4.,    3.,    0.,    0.,   -6./ 1940
      data g9/ -30594.,-2285., 5810.,-1244., 2990.,-1702., 1578.,  477., 1945
     1           1282.,-1834., -499., 1255.,  186.,  913.,  -11.,  944., 1945
     2            776.,  144.,  544., -276., -421.,  -55.,  304., -178., 1945
     3           -253.,  346.,  -12.,  194.,   95.,  -20.,  -67., -142., 1945
     4           -119.,  -82.,   82.,   59.,   57.,    6.,    6.,  100., 1945
     5           -246.,   16.,  -25.,   -9.,   21.,  -16., -104.,  -39., 1945
     6             70.,  -40.,  -45.,    0.,  -18.,    0.,    2.,  -29., 1945
     7              6.,  -10.,   28.,   15.,  -17.,   29.,  -22.,   13., 1945
     8              7.,   12.,   -8.,  -21.,   -5.,  -12.,    9.,   -7., 1945
     9              7.,    2.,  -10.,   18.,    7.,    3.,    2.,  -11., 1945
     a              5.,  -21.,  -27.,    1.,   17.,  -11.,   29.,    3., 1945
     b             -9.,   16.,    4.,   -3.,    9.,   -4.,    6.,   -3., 1945
     c              1.,   -4.,    8.,   -3.,   11.,    5.,    1.,    1., 1945
     d              2.,  -20.,   -5.,   -1.,   -1.,   -6.,    8.,    6., 1945
     e             -1.,   -4.,   -3.,   -2.,    5.,    0.,   -2.,   -2./ 1945
      data ga/ -30554.,-2250., 5815.,-1341., 2998.,-1810., 1576.,  381., 1950
     1           1297.,-1889., -476., 1274.,  206.,  896.,  -46.,  954., 1950
     2            792.,  136.,  528., -278., -408.,  -37.,  303., -210., 1950
     3           -240.,  349.,    3.,  211.,  103.,  -20.,  -87., -147., 1950
     4           -122.,  -76.,   80.,   54.,   57.,   -1.,    4.,   99., 1950
     5           -247.,   33.,  -16.,  -12.,   12.,  -12., -105.,  -30., 1950
     6             65.,  -55.,  -35.,    2.,  -17.,    1.,    0.,  -40., 1950
     7             10.,   -7.,   36.,    5.,  -18.,   19.,  -16.,   22., 1950
     8             15.,    5.,   -4.,  -22.,   -1.,    0.,   11.,  -21., 1950
     9             15.,   -8.,  -13.,   17.,    5.,   -4.,   -1.,  -17., 1950
     a              3.,   -7.,  -24.,   -1.,   19.,  -25.,   12.,   10., 1950
     b              2.,    5.,    2.,   -5.,    8.,   -2.,    8.,    3., 1950
     c            -11.,    8.,   -7.,   -8.,    4.,   13.,   -1.,   -2., 1950
     d             13.,  -10.,   -4.,    2.,    4.,   -3.,   12.,    6., 1950
     e              3.,   -3.,    2.,    6.,   10.,   11.,    3.,    8./ 1950
      data gb/ -30500.,-2215., 5820.,-1440., 3003.,-1898., 1581.,  291., 1955
     1           1302.,-1944., -462., 1288.,  216.,  882.,  -83.,  958., 1955
     2            796.,  133.,  510., -274., -397.,  -23.,  290., -230., 1955
     3           -229.,  360.,   15.,  230.,  110.,  -23.,  -98., -152., 1955
     4           -121.,  -69.,   78.,   47.,   57.,   -9.,    3.,   96., 1955
     5           -247.,   48.,   -8.,  -16.,    7.,  -12., -107.,  -24., 1955
     6             65.,  -56.,  -50.,    2.,  -24.,   10.,   -4.,  -32., 1955
     7              8.,  -11.,   28.,    9.,  -20.,   18.,  -18.,   11., 1955
     8              9.,   10.,   -6.,  -15.,  -14.,    5.,    6.,  -23., 1955
     9             10.,    3.,   -7.,   23.,    6.,   -4.,    9.,  -13., 1955
     a              4.,    9.,  -11.,   -4.,   12.,   -5.,    7.,    2., 1955
     b              6.,    4.,   -2.,    1.,   10.,    2.,    7.,    2., 1955
     c             -6.,    5.,    5.,   -3.,   -5.,   -4.,   -1.,    0., 1955
     d              2.,   -8.,   -3.,   -2.,    7.,   -4.,    4.,    1., 1955
     e             -2.,   -3.,    6.,    7.,   -2.,   -1.,    0.,   -3./ 1955
      data gc/ -30421.,-2169., 5791.,-1555., 3002.,-1967., 1590.,  206., 1960
     1           1302.,-1992., -414., 1289.,  224.,  878., -130.,  957., 1960
     2            800.,  135.,  504., -278., -394.,    3.,  269., -255., 1960
     3           -222.,  362.,   16.,  242.,  125.,  -26., -117., -156., 1960
     4           -114.,  -63.,   81.,   46.,   58.,  -10.,    1.,   99., 1960
     5           -237.,   60.,   -1.,  -20.,   -2.,  -11., -113.,  -17., 1960
     6             67.,  -56.,  -55.,    5.,  -28.,   15.,   -6.,  -32., 1960
     7              7.,   -7.,   23.,   17.,  -18.,    8.,  -17.,   15., 1960
     8              6.,   11.,   -4.,  -14.,  -11.,    7.,    2.,  -18., 1960
     9             10.,    4.,   -5.,   23.,   10.,    1.,    8.,  -20., 1960
     a              4.,    6.,  -18.,    0.,   12.,   -9.,    2.,    1., 1960
     b              0.,    4.,   -3.,   -1.,    9.,   -2.,    8.,    3., 1960
     c              0.,   -1.,    5.,    1.,   -3.,    4.,    4.,    1., 1960
     d              0.,    0.,   -1.,    2.,    4.,   -5.,    6.,    1., 1960
     e              1.,   -1.,   -1.,    6.,    2.,    0.,    0.,   -7./ 1960
      data gd/ -30334.,-2119., 5776.,-1662., 2997.,-2016., 1594.,  114., 1965
     1           1297.,-2038., -404., 1292.,  240.,  856., -165.,  957., 1965
     2            804.,  148.,  479., -269., -390.,   13.,  252., -269., 1965
     3           -219.,  358.,   19.,  254.,  128.,  -31., -126., -157., 1965
     4            -97.,  -62.,   81.,   45.,   61.,  -11.,    8.,  100., 1965
     5           -228.,   68.,    4.,  -32.,    1.,   -8., -111.,   -7., 1965
     6             75.,  -57.,  -61.,    4.,  -27.,   13.,   -2.,  -26., 1965
     7              6.,   -6.,   26.,   13.,  -23.,    1.,  -12.,   13., 1965
     8              5.,    7.,   -4.,  -12.,  -14.,    9.,    0.,  -16., 1965
     9              8.,    4.,   -1.,   24.,   11.,   -3.,    4.,  -17., 1965
     a              8.,   10.,  -22.,    2.,   15.,  -13.,    7.,   10., 1965
     b             -4.,   -1.,   -5.,   -1.,   10.,    5.,   10.,    1., 1965
     c             -4.,   -2.,    1.,   -2.,   -3.,    2.,    2.,    1., 1965
     d             -5.,    2.,   -2.,    6.,    4.,   -4.,    4.,    0., 1965
     e              0.,   -2.,    2.,    3.,    2.,    0.,    0.,   -6./ 1965
      data ge/ -30220.,-2068., 5737.,-1781., 3000.,-2047., 1611.,   25., 1970
     1           1287.,-2091., -366., 1278.,  251.,  838., -196.,  952., 1970
     2            800.,  167.,  461., -266., -395.,   26.,  234., -279., 1970
     3           -216.,  359.,   26.,  262.,  139.,  -42., -139., -160., 1970
     4            -91.,  -56.,   83.,   43.,   64.,  -12.,   15.,  100., 1970
     5           -212.,   72.,    2.,  -37.,    3.,   -6., -112.,    1., 1970
     6             72.,  -57.,  -70.,    1.,  -27.,   14.,   -4.,  -22., 1970
     7              8.,   -2.,   23.,   13.,  -23.,   -2.,  -11.,   14., 1970
     8              6.,    7.,   -2.,  -15.,  -13.,    6.,   -3.,  -17., 1970
     9              5.,    6.,    0.,   21.,   11.,   -6.,    3.,  -16., 1970
     a              8.,   10.,  -21.,    2.,   16.,  -12.,    6.,   10., 1970
     b             -4.,   -1.,   -5.,    0.,   10.,    3.,   11.,    1., 1970
     c             -2.,   -1.,    1.,   -3.,   -3.,    1.,    2.,    1., 1970
     d             -5.,    3.,   -1.,    4.,    6.,   -4.,    4.,    0., 1970
     e              1.,   -1.,    0.,    3.,    3.,    1.,   -1.,   -4./ 1970
      data gf/ -30100.,-2013., 5675.,-1902., 3010.,-2067., 1632.,  -68., 1975
     1           1276.,-2144., -333., 1260.,  262.,  830., -223.,  946., 1975
     2            791.,  191.,  438., -265., -405.,   39.,  216., -288., 1975
     3           -218.,  356.,   31.,  264.,  148.,  -59., -152., -159., 1975
     4            -83.,  -49.,   88.,   45.,   66.,  -13.,   28.,   99., 1975
     5           -198.,   75.,    1.,  -41.,    6.,   -4., -111.,   11., 1975
     6             71.,  -56.,  -77.,    1.,  -26.,   16.,   -5.,  -14., 1975
     7             10.,    0.,   22.,   12.,  -23.,   -5.,  -12.,   14., 1975
     8              6.,    6.,   -1.,  -16.,  -12.,    4.,   -8.,  -19., 1975
     9              4.,    6.,    0.,   18.,   10.,  -10.,    1.,  -17., 1975
     a              7.,   10.,  -21.,    2.,   16.,  -12.,    7.,   10., 1975
     b             -4.,   -1.,   -5.,   -1.,   10.,    4.,   11.,    1., 1975
     c             -3.,   -2.,    1.,   -3.,   -3.,    1.,    2.,    1., 1975
     d             -5.,    3.,   -2.,    4.,    5.,   -4.,    4.,   -1., 1975
     e              1.,   -1.,    0.,    3.,    3.,    1.,   -1.,   -5./ 1975
      data gg/ -29992.,-1956., 5604.,-1997., 3027.,-2129., 1663., -200., 1980
     1           1281.,-2180., -336., 1251.,  271.,  833., -252.,  938., 1980
     2            782.,  212.,  398., -257., -419.,   53.,  199., -297., 1980
     3           -218.,  357.,   46.,  261.,  150.,  -74., -151., -162., 1980
     4            -78.,  -48.,   92.,   48.,   66.,  -15.,   42.,   93., 1980
     5           -192.,   71.,    4.,  -43.,   14.,   -2., -108.,   17., 1980
     6             72.,  -59.,  -82.,    2.,  -27.,   21.,   -5.,  -12., 1980
     7             16.,    1.,   18.,   11.,  -23.,   -2.,  -10.,   18., 1980
     8              6.,    7.,    0.,  -18.,  -11.,    4.,   -7.,  -22., 1980
     9              4.,    9.,    3.,   16.,    6.,  -13.,   -1.,  -15., 1980
     a              5.,   10.,  -21.,    1.,   16.,  -12.,    9.,    9., 1980
     b             -5.,   -3.,   -6.,   -1.,    9.,    7.,   10.,    2., 1980
     c             -6.,   -5.,    2.,   -4.,   -4.,    1.,    2.,    0., 1980
     d             -5.,    3.,   -2.,    6.,    5.,   -4.,    3.,    0., 1980
     e              1.,   -1.,    2.,    4.,    3.,    0.,    0.,   -6./ 1980
      data gi/ -29873.,-1905., 5500.,-2072., 3044.,-2197., 1687., -306., 1985
     1           1296.,-2208., -310., 1247.,  284.,  829., -297.,  936., 1985
     2            780.,  232.,  361., -249., -424.,   69.,  170., -297., 1985
     3           -214.,  355.,   47.,  253.,  150.,  -93., -154., -164., 1985
     4            -75.,  -46.,   95.,   53.,   65.,  -16.,   51.,   88., 1985
     5           -185.,   69.,    4.,  -48.,   16.,   -1., -102.,   21., 1985
     6             74.,  -62.,  -83.,    3.,  -27.,   24.,   -2.,   -6., 1985
     7             20.,    4.,   17.,   10.,  -23.,    0.,   -7.,   21., 1985
     8              6.,    8.,    0.,  -19.,  -11.,    5.,   -9.,  -23., 1985
     9              4.,   11.,    4.,   14.,    4.,  -15.,   -4.,  -11., 1985
     a              5.,   10.,  -21.,    1.,   15.,  -12.,    9.,    9., 1985
     b             -6.,   -3.,   -6.,   -1.,    9.,    7.,    9.,    1., 1985
     c             -7.,   -5.,    2.,   -4.,   -4.,    1.,    3.,    0., 1985
     d             -5.,    3.,   -2.,    6.,    5.,   -4.,    3.,    0., 1985
     e              1.,   -1.,    2.,    4.,    3.,    0.,    0.,   -6./ 1985
      data gj/ -29775.,-1848., 5406.,-2131., 3059.,-2279., 1686., -373., 1990
     1           1314.,-2239., -284., 1248.,  293.,  802., -352.,  939., 1990
     2            780.,  247.,  325., -240., -423.,   84.,  141., -299., 1990
     3           -214.,  353.,   46.,  245.,  154., -109., -153., -165., 1990
     4            -69.,  -36.,   97.,   61.,   65.,  -16.,   59.,   82., 1990
     5           -178.,   69.,    3.,  -52.,   18.,    1.,  -96.,   24., 1990
     6             77.,  -64.,  -80.,    2.,  -26.,   26.,    0.,   -1., 1990
     7             21.,    5.,   17.,    9.,  -23.,    0.,   -4.,   23., 1990
     8              5.,   10.,   -1.,  -19.,  -10.,    6.,  -12.,  -22., 1990
     9              3.,   12.,    4.,   12.,    2.,  -16.,   -6.,  -10., 1990
     a              4.,    9.,  -20.,    1.,   15.,  -12.,   11.,    9., 1990
     b             -7.,   -4.,   -7.,   -2.,    9.,    7.,    8.,    1., 1990
     c             -7.,   -6.,    2.,   -3.,   -4.,    2.,    2.,    1., 1990
     d             -5.,    3.,   -2.,    6.,    4.,   -4.,    3.,    0., 1990
     e              1.,   -2.,    3.,    3.,    3.,   -1.,    0.,   -6./ 1990
      data gk/ -29692.,-1784., 5306.,-2200., 3070.,-2366., 1681., -413., 1995
     1           1335.,-2267., -262., 1249.,  302.,  759., -427.,  940., 1995
     2            780.,  262.,  290., -236., -418.,   97.,  122., -306., 1995
     3           -214.,  352.,   46.,  235.,  165., -118., -143., -166., 1995
     4            -55.,  -17.,  107.,   68.,   67.,  -17.,   68.,   72., 1995
     5           -170.,   67.,   -1.,  -58.,   19.,    1.,  -93.,   36., 1995
     6             77.,  -72.,  -69.,    1.,  -25.,   28.,    4.,    5., 1995
     7             24.,    4.,   17.,    8.,  -24.,   -2.,   -6.,   25., 1995
     8              6.,   11.,   -6.,  -21.,   -9.,    8.,  -14.,  -23., 1995
     9              9.,   15.,    6.,   11.,   -5.,  -16.,   -7.,   -4., 1995
     a              4.,    9.,  -20.,    3.,   15.,  -10.,   12.,    8., 1995
     b             -6.,   -8.,   -8.,   -1.,    8.,   10.,    5.,   -2., 1995
     c             -8.,   -8.,    3.,   -3.,   -6.,    1.,    2.,    0., 1995
     d             -4.,    4.,   -1.,    5.,    4.,   -5.,    2.,   -1., 1995
     e              2.,   -2.,    5.,    1.,    1.,   -2.,    0.,   -7., 1995
     f           75*0./                                                  1995
      data gl/ -29619.4,-1728.2, 5186.1,-2267.7, 3068.4,-2481.6, 1670.9, 2000
     1           -458.0, 1339.6,-2288.0, -227.6, 1252.1,  293.4,  714.5, 2000
     2           -491.1,  932.3,  786.8,  272.6,  250.0, -231.9, -403.0, 2000
     3            119.8,  111.3, -303.8, -218.8,  351.4,   43.8,  222.3, 2000
     4            171.9, -130.4, -133.1, -168.6,  -39.3,  -12.9,  106.3, 2000
     5             72.3,   68.2,  -17.4,   74.2,   63.7, -160.9,   65.1, 2000
     6             -5.9,  -61.2,   16.9,    0.7,  -90.4,   43.8,   79.0, 2000
     7            -74.0,  -64.6,    0.0,  -24.2,   33.3,    6.2,    9.1, 2000
     8             24.0,    6.9,   14.8,    7.3,  -25.4,   -1.2,   -5.8, 2000
     9             24.4,    6.6,   11.9,   -9.2,  -21.5,   -7.9,    8.5, 2000
     a            -16.6,  -21.5,    9.1,   15.5,    7.0,    8.9,   -7.9, 2000
     b            -14.9,   -7.0,   -2.1,    5.0,    9.4,  -19.7,    3.0, 2000
     c             13.4,   -8.4,   12.5,    6.3,   -6.2,   -8.9,   -8.4, 2000
     d             -1.5,    8.4,    9.3,    3.8,   -4.3,   -8.2,   -8.2, 2000
     e              4.8,   -2.6,   -6.0,    1.7,    1.7,    0.0,   -3.1, 2000
     f              4.0,   -0.5,    4.9,    3.7,   -5.9,    1.0,   -1.2, 2000
     g              2.0,   -2.9,    4.2,    0.2,    0.3,   -2.2,   -1.1, 2000
     h             -7.4,    2.7,   -1.7,    0.1,   -1.9,    1.3,    1.5, 2000
     i             -0.9,   -0.1,   -2.6,    0.1,    0.9,   -0.7,   -0.7, 2000
     j              0.7,   -2.8,    1.7,   -0.9,    0.1,   -1.2,    1.2, 2000
     k             -1.9,    4.0,   -0.9,   -2.2,   -0.3,   -0.4,    0.2, 2000
     l              0.3,    0.9,    2.5,   -0.2,   -2.6,    0.9,    0.7, 2000
     m             -0.5,    0.3,    0.3,    0.0,   -0.3,    0.0,   -0.4, 2000
     n              0.3,   -0.1,   -0.9,   -0.2,   -0.4,   -0.4,    0.8, 2000
     o             -0.2,   -0.9,   -0.9,    0.3,    0.2,    0.1,    1.8, 2000
     p             -0.4,   -0.4,    1.3,   -1.0,   -0.4,   -0.1,    0.7, 2000
     q              0.7,   -0.4,    0.3,    0.3,    0.6,   -0.1,    0.3, 2000
     r              0.4,   -0.2,    0.0,   -0.5,    0.1,   -0.9/         2000
      data gm/-29554.63,-1669.05, 5077.99,-2337.24, 3047.69,-2594.50,    2005
     1          1657.76, -515.43, 1336.30,-2305.83, -198.86, 1246.39,    2005
     2           269.72,  672.51, -524.72,  920.55,  797.96,  282.07,    2005
     3           210.65, -225.23, -379.86,  145.15,  100.00, -305.36,    2005
     4          -227.00,  354.41,   42.72,  208.95,  180.25, -136.54,    2005
     5          -123.45, -168.05,  -19.57,  -13.55,  103.85,   73.60,    2005
     6            69.56,  -20.33,   76.74,   54.75, -151.34,   63.63,    2005
     7           -14.58,  -63.53,   14.58,    0.24,  -86.36,   50.94,    2005
     8            79.88,  -74.46,  -61.14,   -1.65,  -22.57,   38.73,    2005
     9             6.82,   12.30,   25.35,    9.37,   10.93,    5.42,    2005
     a           -26.32,    1.94,   -4.64,   24.80,    7.62,   11.20,    2005
     b           -11.73,  -20.88,   -6.88,    9.83,  -18.11,  -19.71,    2005
     c            10.17,   16.22,    9.36,    7.61,  -11.25,  -12.76,    2005
     d            -4.87,   -0.06,    5.58,    9.76,  -20.11,    3.58,    2005
     e            12.69,   -6.94,   12.67,    5.01,   -6.72,  -10.76,    2005
     f            -8.16,   -1.25,    8.10,    8.76,    2.92,   -6.66,    2005
     g            -7.73,   -9.22,    6.01,   -2.17,   -6.12,    2.19,    2005
     h             1.42,    0.10,   -2.35,    4.46,   -0.15,    4.76,    2005
     i             3.06,   -6.58,    0.29,   -1.01,    2.06,   -3.47,    2005
     j             3.77,   -0.86,   -0.21,   -2.31,   -2.09,   -7.93,    2005
     k             2.95,   -1.60,    0.26,   -1.88,    1.44,    1.44,    2005
     l            -0.77,   -0.31,   -2.27,    0.29,    0.90,   -0.79,    2005
     m            -0.58,    0.53,   -2.69,    1.80,   -1.08,    0.16,    2005
     n            -1.58,    0.96,   -1.90,    3.99,   -1.39,   -2.15,    2005
     o            -0.29,   -0.55,    0.21,    0.23,    0.89,    2.38,    2005
     p            -0.38,   -2.63,    0.96,    0.61,   -0.30,    0.40,    2005
     q             0.46,    0.01,   -0.35,    0.02,   -0.36,    0.28,    2005
     r             0.08,   -0.87,   -0.49,   -0.34,   -0.08,    0.88,    2005
     s            -0.16,   -0.88,   -0.76,    0.30,    0.33,    0.28,    2005
     t             1.72,   -0.43,   -0.54,    1.18,   -1.07,   -0.37,    2005
     u            -0.04,    0.75,    0.63,   -0.26,    0.21,    0.35,    2005
     v             0.53,   -0.05,    0.38,    0.41,   -0.22,   -0.10,    2005
     w            -0.57,   -0.18,   -0.82/                               2005
      data gp/-29496.5,-1585.9, 4945.1,-2396.6, 3026.0,-2707.7, 1668.6,  2010
     1          -575.4, 1339.7,-2326.3, -160.5, 1231.7,  251.7,  634.2,  2010
     2          -536.8,  912.6,  809.0,  286.4,  166.6, -211.2, -357.1,  2010
     3           164.4,   89.7, -309.2, -231.1,  357.2,   44.7,  200.3,  2010
     4           188.9, -141.2, -118.1, -163.1,    0.1,   -7.7,  100.9,  2010
     5            72.8,   68.6,  -20.8,   76.0,   44.2, -141.4,   61.5,  2010
     6           -22.9,  -66.3,   13.1,    3.1,  -77.9,   54.9,   80.4,  2010
     7           -75.0,  -57.8,   -4.7,  -21.2,   45.3,    6.6,   14.0,  2010
     8            24.9,   10.4,    7.0,    1.6,  -27.7,    4.9,   -3.4,  2010
     9            24.3,    8.2,   10.9,  -14.5,  -20.0,   -5.7,   11.9,  2010
     a           -19.3,  -17.4,   11.6,   16.7,   10.9,    7.1,  -14.1,  2010
     b           -10.8,   -3.7,    1.7,    5.4,    9.4,  -20.5,    3.4,  2010
     c            11.6,   -5.3,   12.8,    3.1,   -7.2,  -12.4,   -7.4,  2010
     d            -0.8,    8.0,    8.4,    2.2,   -8.4,   -6.1,  -10.1,  2010
     e             7.0,   -2.0,   -6.3,    2.8,    0.9,   -0.1,   -1.1,  2010
     f             4.7,   -0.2,    4.4,    2.5,   -7.2,   -0.3,   -1.0,  2010
     g             2.2,   -4.0,    3.1,   -2.0,   -1.0,   -2.0,   -2.8,  2010
     h            -8.3,    3.0,   -1.5,    0.1,   -2.1,    1.7,    1.6,  2010
     i            -0.6,   -0.5,   -1.8,    0.5,    0.9,   -0.8,   -0.4,  2010
     j             0.4,   -2.5,    1.8,   -1.3,    0.2,   -2.1,    0.8,  2010
     k            -1.9,    3.8,   -1.8,   -2.1,   -0.2,   -0.8,    0.3,  2010
     l             0.3,    1.0,    2.2,   -0.7,   -2.5,    0.9,    0.5,  2010
     m            -0.1,    0.6,    0.5,    0.0,   -0.4,    0.1,   -0.4,  2010
     n             0.3,    0.2,   -0.9,   -0.8,   -0.2,    0.0,    0.8,  2010
     o            -0.2,   -0.9,   -0.8,    0.3,    0.3,    0.4,    1.7,  2010
     p            -0.4,   -0.6,    1.1,   -1.2,   -0.3,   -0.1,    0.8,  2010
     q             0.5,   -0.2,    0.1,    0.4,    0.5,    0.0,    0.4,  2010
     r             0.4,   -0.2,   -0.3,   -0.5,   -0.3,   -0.8/          2010
      data gq/    11.4,   16.7,  -28.8,  -11.3,   -3.9,  -23.0,    2.7,  2012
     1           -12.9,    1.3,   -3.9,    8.6,   -2.9,   -2.9,   -8.1,  2012
     2            -2.1,   -1.4,    2.0,    0.4,   -8.9,    3.2,    4.4,  2012
     3             3.6,   -2.3,   -0.8,   -0.5,    0.5,    0.5,   -1.5,  2012
     4             1.5,   -0.7,    0.9,    1.3,    3.7,    1.4,   -0.6,  2012
     5            -0.3,   -0.3,   -0.1,   -0.3,   -2.1,    1.9,   -0.4,  2012
     6            -1.6,   -0.5,   -0.2,    0.8,    1.8,    0.5,    0.2,  2012
     7            -0.1,    0.6,   -0.6,    0.3,    1.4,   -0.2,    0.3,  2012
     8            -0.1,    0.1,   -0.8,   -0.8,   -0.3,    0.4,    0.2,  2012
     9            -0.1,    0.1,    0.0,   -0.5,    0.2,    0.3,    0.5,  2012
     a            -0.3,    0.4,    0.3,    0.1,    0.2,   -0.1,   -0.5,  2012
     b             0.4,    0.2,    0.4,115*0.0/                          2012
c
c     set initial values
c
      x     = 0.0
      y     = 0.0
      z     = 0.0
      if (date.lt.1900.0.or.date.gt.2020.0) go to 11
      if (date.gt.2015.0) write (6,960) date
  960 format (/' This version of the IGRF is intended for use up',
     1        ' to 2015.0.'/' values for',f9.3,' will be computed',
     2        ' but may be of reduced accuracy'/)
      if (date.ge.2010.0) go to 1
      t     = 0.2*(date - 1900.0)
      ll    = t
      one   = ll
      t     = t - one
c
c     SH models before 1995.0 are only to degree 10
c
      if (date.lt.1995.0) then
       nmx   = 10
       nc    = nmx*(nmx+2)
       ll    = nc*ll
       kmx   = (nmx+1)*(nmx+2)/2
      else
       nmx   = 13
       nc    = nmx*(nmx+2)
       ll    = 0.2*(date - 1995.0)
c
c     19 is the number of SH models that extend to degree 10
c
       ll    = 120*19 + nc*ll
       kmx   = (nmx+1)*(nmx+2)/2
      endif
      tc    = 1.0 - t
      if (isv.eq.1) then
       tc = -0.2
       t = 0.2
      end if
      go to 2
c
    1 t     = date - 2010.0
      tc    = 1.0
      if (isv.eq.1) then
       t = 1.0
       tc = 0.0
      end if
c
c     pointer for last coefficient in pen-ultimate set of MF coefficients...
c
      ll    = 2865
      nmx   = 13
      nc    = nmx*(nmx+2)
      kmx   = (nmx+1)*(nmx+2)/2
    2 r     = alt
      one   = colat*0.017453292
      ct    = cos(one)
      st    = sin(one)
      one   = elong*0.017453292
      cl(1) = cos(one)
      sl(1) = sin(one)
      cd    = 1.0
      sd    = 0.0
      l     = 1
      m     = 1
      n     = 0
      if (itype.eq.2) go to 3
c
c     conversion from geodetic to geocentric coordinates
c     (using the WGS84 spheroid)
c
      a2    = 40680631.6
      b2    = 40408296.0
      one   = a2*st*st
      two   = b2*ct*ct
      three = one + two
      rho   = sqrt(three)
      r     = sqrt(alt*(alt + 2.0*rho) + (a2*one + b2*two)/three)
      cd    = (alt + rho)/r
      sd    = (a2 - b2)/rho*ct*st/r
      one   = ct
      ct    = ct*cd -  st*sd
      st    = st*cd + one*sd
c
    3 ratio = 6371.2/r
      rr    = ratio*ratio
c
c     computation of Schmidt quasi-normal coefficients p and x(=q)
c
      p(1)  = 1.0
      p(3)  = st
      q(1)  = 0.0
      q(3)  =  ct
      do 10 k=2,kmx
       if (n.ge.m) go to 4
       m     = 0
       n     = n + 1
       rr    = rr*ratio
       fn    = n
       gn    = n - 1
    4  fm    = m
       if (m.ne.n) go to 5
       if (k.eq.3) go to 6
       one   = sqrt(1.0 - 0.5/fm)
       j     = k - n - 1
       p(k)  = one*st*p(j)
       q(k)  = one*(st*q(j) + ct*p(j))
       cl(m) = cl(m-1)*cl(1) - sl(m-1)*sl(1)
       sl(m) = sl(m-1)*cl(1) + cl(m-1)*sl(1)
       go to 6
    5  gmm    = m*m
       one   = sqrt(fn*fn - gmm)
       two   = sqrt(gn*gn - gmm)/one
       three = (fn + gn)/one
       i     = k - n
       j     = i - n + 1
       p(k)  = three*ct*p(i) - two*p(j)
       q(k)  = three*(ct*q(i) - st*p(i)) - two*q(j)
c
c     synthesis of x, y and z in geocentric coordinates
c
    6  lm    = ll + l
       one   = (tc*gh(lm) + t*gh(lm+nc))*rr
       if (m.eq.0) go to 9
       two   = (tc*gh(lm+1) + t*gh(lm+nc+1))*rr
       three = one*cl(m) + two*sl(m)
       x     = x + three*q(k)
       z     = z - (fn + 1.0)*three*p(k)
       if (st.eq.0.0) go to 7
       y     = y + (one*sl(m) - two*cl(m))*fm*p(k)/st
       go to 8
    7  y     = y + (one*sl(m) - two*cl(m))*q(k)*ct
    8  l     = l + 2
       go to 10
    9  x     = x + one*q(k)
       z     = z - (fn + 1.0)*one*p(k)
       l     = l + 1
   10 m     = m + 1
c
c     conversion to coordinate system specified by itype
c
      one   = x
      x     = x*cd +   z*sd
      z     = z*cd - one*sd
      f     = sqrt(x*x + y*y + z*z)
c
      return
c
c     error return if date out of bounds
c
   11 f     = 1.0d8
      write (6,961) date
  961 format (/' This subroutine will not work with a date of',
     1        f9.3,'.  Date must be in the range 1900.0.ge.date',
     2        '.le.2020.0. On return f = 1.0d8., x = y = z = 0.')
      return
      end

k	end module

max_stars_repo_name: DearDon/Geomagnetics

max_issues_repo_name: DearDon/Geomagnetics

max_forks_repo_name: DearDon/Geomagnetics


text:
   PROGRAM GALAHAD_TRS_EXAMPLE  !  GALAHAD 2.2 - 05/06/2008 AT 13:30 GMT.
   USE GALAHAD_TRS_DOUBLE                          ! double precision version
   IMPLICIT NONE
   INTEGER, PARAMETER :: wp = KIND( 1.0D+0 )       ! set precision
   REAL ( KIND = wp ), PARAMETER :: one = 1.0_wp
   INTEGER, PARAMETER :: n = 3, h_ne = 4           ! problem dimensions
   INTEGER :: s
   REAL ( KIND = wp ), DIMENSION( n ) :: C, X
   TYPE ( SMT_type ) :: H
   TYPE ( TRS_data_type ) :: data
   TYPE ( TRS_control_type ) :: control
   TYPE ( TRS_inform_type ) :: inform
   REAL ( KIND = wp ) :: f = one                ! constant term, f
   REAL ( KIND = wp ) :: radius = one           ! trust-region radius
   C = (/ 0.0_wp, 2.0_wp, 0.0_wp /)             ! linear term, c
!  C = (/ 0.0_wp, 2.0_wp, 0.0001_wp /)             ! linear term, c
!  C = (/ 5.0_wp, 0.0_wp, 4.0_wp /)             ! linear term, c
   CALL SMT_put( H%type, 'COORDINATE', s )      ! Specify co-ordinate
   ALLOCATE( H%val( h_ne ), H%row( h_ne ), H%col( h_ne ) )
   H%val = (/ 1.0_wp, 2.0_wp, 3.0_wp, 4.0_wp /) ! Hessian, H
!  H%val = (/ 1.0_wp, 2.0_wp, 3.0_wp, 1.0_wp /) ! Hessian, H
   H%row = (/ 1, 2, 3, 3 /)                     ! NB lower triangle
   H%col = (/ 1, 2, 3, 1 /) ; H%ne = h_ne
   CALL TRS_initialize( data, control )         ! Initialize control parameters
!  control%print_level = 1
!  control%initial_multiplier = 4.0_wp
   CALL TRS_solve( n, radius, f, C, H, X, data, control, inform ) ! solve problem
   IF ( inform%status == 0 ) THEN !  Successful return
     WRITE( 6, "( 1X, I0, ' factorizations. Solution and Lagrange multiplier =',&
    &    2ES12.4 )" ) inform%factorizations, inform%obj, inform%multiplier
   ELSE  !  Error returns
     WRITE( 6, "( ' TRS_solve exit status = ', I0 ) " ) inform%status
   END IF
   CALL TRS_terminate( data, control, inform )  ! delete internal workspace
   END PROGRAM GALAHAD_TRS_EXAMPLE

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
C (C) Copyright 1996-2016 ECMWF.
C
C This software is licensed under the terms of the Apache Licence Version 2.0
C which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
C In applying this licence, ECMWF does not waive the privileges and immunities
C granted to it by virtue of its status as an intergovernmental organisation nor
C does it submit to any jurisdiction.
C

      program polys

      call popen
c      call psetc('legend_box_mode', 'positional')
c      call psetr('legend_box_x_position', 25.)
c      call psetr('legend_box_y_position', 7.)
c      call psetr('legend_box_x_length', 5.)
c      call psetr('legend_box_y_length', 10.)
c     call psetc('legend_box_blanking', 'on')
c      call psetc('legend_box_border', 'on')
      call psetc('legend_title', 'on')
      call psetc('legend_title_text', 'My title')
      call psetc('legend_display_type', 'continuous')
      call psetc('legend_text_format', '(f5.2)')

        call psetr ('subpage_lower_left_longitude',  -30.0)
        call psetr ('subpage_upper_right_longitude',  77.0)
        call psetr ('subpage_upper_right_latitude',   75.0)
        call psetr ('subpage_lower_left_latitude',    25.0)
        call psetc('legend', 'on')

        call psetc('polyline_input_positions_filename',
     +         '../data/catchment.pos')
        call psetc('polyline_input_values_filename',
     +         '../data/catchment.val')
        call psetr('polyline_input_break_indicator', -999.)
        call psetc('polyline_shade', 'on')
        call psetc('polyline_shade_colour_method',  'calculate')
        call psetc('polyline_shade_max_level_colour', 'red')

        call psetc('polyline_shade_min_level_colour', 'blue')
        call psetc('polyline_shade_colour_direction', 'clockwise')
        call pline
      call pcoast
       call pnew("page")
      call psetc('legend_box_mode', 'positional')
      call psetr('legend_box_x_position', 25.)
      call psetr('legend_box_y_position', 7.)
      call psetr('legend_box_x_length', 5.)
       call psetr('polyline_shade_max_level', 8.)
        call psetr('polyline_shade_min_level', 4.)
      call psetr('legend_box_y_length', 10.)
      call psetc('legend_box_blanking', 'on')
      call psetc('legend_box_border', 'on')
       call psetc('polyline_shade_level_selection_type', 'interval')
       call psetr('polyline_interval', 0.5)

      call pline
      call pcoast

      call pclose

      end

max_stars_repo_name: dtip/magics

max_issues_repo_name: dtip/magics

max_forks_repo_name: dtip/magics


text:
      subroutine cpu_second (cpu,user,sys)
      real cpu,user,sys
c
c  ------------------------------------------------
c  returns elapsed cpu time since start of job (sec)
c  ------------------------------------------------
c
c nil version; returns -1 for all clocks; can be used when there is no way
c to measure CPU time
c
      cpu = -1.
      user = -1.
      sys = -1.
c
      return
      end

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
MODULE fitChi2_external

  USE auxil, ONLY: parinfo_type, indpar_type, Qabs_type
  USE utilities, ONLY: DP
  IMPLICIT NONE
  PRIVATE

  INTEGER, SAVE, PUBLIC :: Nx, Ny, NwOBS, xOBS, yOBS, jw
  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE, PUBLIC        :: iwfree
  REAL(DP), DIMENSION(:), ALLOCATABLE, SAVE, PUBLIC       :: wOBS, nuOBS
  REAL(DP), DIMENSION(:,:), ALLOCATABLE, SAVE, PUBLIC :: extinct
  REAL(DP), DIMENSION(:), ALLOCATABLE, SAVE, PUBLIC       :: Fnu_mod, resid
  REAL(DP), DIMENSION(:,:,:), ALLOCATABLE, SAVE, PUBLIC   :: FnuOBS, dFnuOBS
  REAL(DP), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE, PUBLIC :: invLcovarOBS
  LOGICAL, SAVE, PUBLIC :: iid ! indpdt identically dist. (wvl)
  LOGICAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE, PUBLIC    :: mask

  !! Init param
  !!------------
  TYPE(indpar_type), SAVE, PUBLIC                             :: ind
  TYPE(parinfo_type), DIMENSION(:), ALLOCATABLE, SAVE, PUBLIC :: parinfo
  TYPE(Qabs_type), DIMENSION(:), ALLOCATABLE, SAVE, PUBLIC    :: Qabs

  PUBLIC :: residuals

CONTAINS

  !! Main function for L-M
  !!-----------------------
  FUNCTION residuals(par, NwOBS0)

    USE auxil, ONLY: specModel
    USE utilities, ONLY: DP
    IMPLICIT NONE

    INTEGER, INTENT(IN)                      :: NwOBS0
    REAL(DP), DIMENSION(:), INTENT(IN)       :: par
    INTEGER                                  :: i

    REAL(DP), DIMENSION(NwOBS0)              :: residuals

    Fnu_mod(:) = specModel(wOBS(:), INDPAR=ind, PARVAL=par(:), &
                           QABS=Qabs(:), EXTINCT=extinct(:,:))

    !! Unweighted residuals
    WHERE (mask(xOBS,yOBS,:))
      resid(:) = FnuOBS(xOBS,yOBS,:) - Fnu_mod(:)
    ELSEWHERE
      resid(:) = 0._DP

    END WHERE

    !! Weighted residuals
    ! residuals(:) = MATMUL(invLcovarOBS(xOBS,yOBS,:,:), resid(:))

    !! With keyword iid
    IF (iid) THEN
      FORALL (i=1:NwOBS0) &
        residuals(i) = invLcovarOBS(xOBS,yOBS,i,i) * resid(i)
    ELSE
      residuals(:) = MATMUL(invLcovarOBS(xOBS,yOBS,:,:), resid(:))

    END IF

  END FUNCTION residuals

END MODULE fitChi2_external

!!==========================================================================
!!                           Main Program
!!==========================================================================

PROGRAM test_fitChi2

  USE auxil, ONLY: read_master, specModel, initparam
  USE utilities, ONLY: DP, pring, trimLR, trimeq, timinfo, &
                       banner_program, ustd, isNaN, initiate_clock, time_type
  USE arrays, ONLY: iwhere, closest
  USE constants, ONLY: MKS
  USE statistics, ONLY: MEDIAN
  USE inout, ONLY: read_hdf5, write_hdf5, h5ext, write_ascii, ascext, lenpar, lenpath
  USE grain_optics, ONLY: lendustQ
  USE random, ONLY: generate_newseed
  USE chi2_minimization, ONLY: chi2min_LM
  USE fitChi2_external, ONLY: Nx, Ny, NwOBS, wOBS, nuOBS, FnuOBS, dFnuOBS, &
                              residuals, Fnu_mod, resid, invLcovarOBS, extinct, &
                              ind, parinfo, Qabs, xOBS, yOBS, mask, iwfree, iid
  IMPLICIT NONE

  !! Parameters
  INTEGER, PARAMETER :: ulog = 2
  ! REAL(DP), PARAMETER :: tol = 1.E-10_DP
  CHARACTER(*), PARAMETER :: filOBS = './dat/observations_fitMIR'//h5ext
  CHARACTER(*), PARAMETER :: dirOUT = './out/'
  CHARACTER(*), PARAMETER :: filOUT = dirOUT//'test_fitChi2'//h5ext
  CHARACTER(*), PARAMETER :: fiLOG = './tmp/log_fitChi2'//ascext
  LOGICAL, PARAMETER :: compress = .FALSE.
  LOGICAL, PARAMETER :: debug = .FALSE.

  !! Input variables
  INTEGER :: i, j, x, y, Npar, Nparfree, Nwfree, NiniMC
  INTEGER :: Ncont, Nband, Nline, Nextc, Nstar, Nextra
  CHARACTER(lenpar) :: spec_unit
  CHARACTER(lendustQ), DIMENSION(:), ALLOCATABLE :: labQ
  CHARACTER(lenpar), DIMENSION(:), ALLOCATABLE :: labB, labL
  LOGICAL :: verbose, calib, newseed, newinit, dostop
  INTEGER, DIMENSION(:,:,:), ALLOCATABLE    :: maskint ! convert mask=0 to mask=T

  !! Monte-Carlo parameters
  INTEGER :: ibest
  INTEGER, DIMENSION(:), ALLOCATABLE :: statusMC, NiterMC
  REAL(DP), DIMENSION(:), ALLOCATABLE :: chi2redMC
  REAL(DP), DIMENSION(:,:), ALLOCATABLE :: parerrMC
  REAL(DP), DIMENSION(:,:,:), ALLOCATABLE :: covarMC
  REAL(DP), DIMENSION(:,:,:,:), ALLOCATABLE :: parini

  !! chi2min_LM parameters
  INTEGER, DIMENSION(:,:), ALLOCATABLE      :: status
  INTEGER, DIMENSION(:), ALLOCATABLE        :: itied
  REAL(DP), DIMENSION(:,:), ALLOCATABLE     :: limits!, medSovN
  REAL(DP), DIMENSION(:,:), ALLOCATABLE     :: chi2red
  REAL(DP), DIMENSION(:,:,:), ALLOCATABLE   :: parerr, par
  REAL(DP), DIMENSION(:,:,:,:), ALLOCATABLE :: covarOBS, covpar
  LOGICAL, DIMENSION(:,:), ALLOCATABLE      :: limited, maskS
  ! LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE  :: mask4D

  !! Output variables
  INTEGER, DIMENSION(:,:), ALLOCATABLE :: Niter
  INTEGER, DIMENSION(2) :: unitlog
  TYPE(time_type) :: timestr
  REAL(DP), DIMENSION(:,:,:), ALLOCATABLE :: FnuCONT, &
    FnuBAND, FnuSTAR, Pabs, FnuLINE, FnuMOD
  REAL(DP), DIMENSION(:,:,:,:), ALLOCATABLE :: FnuCONT_tab, &
    FnuBAND_tab, FnuSTAR_tab, Pabs_tab, FnuLINE_tab

  ! Output settings
  !----------------
  CALL INITIATE_CLOCK(timestr)
  OPEN (ulog,FILE=fiLOG,STATUS="REPLACE",ACTION="WRITE")
  unitlog(:) = [ ulog, ustd ]

  !! Wavelengths -> indpdt identically dist. (diagonal covar matrix)
  iid = .TRUE.

  !!------------------------------------------------------------------------
  !!                            I. Read the inputs
  !!------------------------------------------------------------------------

  !! a. Main settings
  !!------------------
  CALL READ_HDF5(wOBS, FILE=filOBS, NAME='Wavelength (microns)', N1=NwOBS)
  CALL READ_MASTER(WAVALL=wOBS(:), &
                   VERBOSE=verbose, NINIMC=NiniMC, &
                   CALIB=calib, NEWSEED=newseed, NEWINIT=newinit, &
                   LABQ=labQ, LABL=labL, LABB=labB, QABS=Qabs, EXTINCT=extinct, &
                   NCONT=Ncont, NBAND=Nband, NLINE=Nline, &
                   NEXTC=Nextc, NSTAR=Nstar, NEXTRA=Nextra, DOSTOP=dostop, &
                   PARINFO=parinfo, INDPAR=ind, NPAR=Npar, SPEC_UNIT=spec_unit)

  IF (newseed) CALL GENERATE_NEWSEED()
  IF (verbose) PRINT*
  DO i=1,MERGE(2,1,verbose)
    CALL BANNER_PROGRAM("LE MIROIR: LEast-squares fitting of Mid-IR emission " &
                        //"OptImized Routine", UNIT=unitlog(i), SWING=.True.)
  END DO

  !! b. Observation sample
  !!-----------------------
  CALL READ_HDF5(DBLARR3D=FnuOBS, FILE=filOBS, &
                 NAME='FnuOBS ('//TRIMLR(spec_unit)//')', N1=Nx, N2=Ny)
  CALL READ_HDF5(DBLARR3D=dFnuOBS, FILE=filOBS, &
                 NAME='dFnuOBS ('//TRIMLR(spec_unit)//')')

  dFnuOBS(:,:,:) = dFnuOBS(:,:,:)! + 0.1*ABS(FnuOBS(:,:,:))
  !! On input, the convention is mask=1 => blocks; mask=0 => passes.
  !! However, within Fortran, the convention is the Fortran mask=T => passes.
  CALL READ_HDF5(INTARR3D=maskint, FILE=filOBS, NAME='NaN mask')
  ! print*, size(maskint, dim=3)
  ALLOCATE(mask(Nx,Ny,NwOBS))
  mask(:,:,:) = ( maskint(:,:,:) == 0 )
  ! print*, count(isnan(dfnuobs) .and. mask)

  ALLOCATE(nuOBS(NwOBS))
  nuOBS(:) = MKS%clight/MKS%micron / wOBS(:)

  !! dfnu positive -> 0 ???
  ! WHERE (isNaN(FnuOBS(:,:,:))) FnuOBS(:,:,:) = 0._DP
  ! WHERE (isNaN(dFnuOBS(:,:,:))) dFnuOBS(:,:,:) = 0._DP

  Nwfree = COUNT( ALL(ALL(.NOT. mask(:,:,:) &
                            .AND. dFnuOBS(:,:,:) > 0._DP,DIM=1),DIM=1) )
  IF (Nwfree > 0) &
    CALL IWHERE(ALL(ALL(.NOT. mask(:,:,:) &
                        .AND. dFnuOBS(:,:,:) > 0._DP,DIM=1),DIM=1), iwfree)

  !! Median signal-to-noise ratio
  ! ALLOCATE(medSovN(Nx,Ny))
  ! medSovN(:,:) = 0._DP
  ! DO x=1,Nx
  !   DO y=1,Ny
  !     IF (ANY(mask(x,y,:))) THEN
  !       medSovN(x,y) = MEDIAN(FnuOBS(x,y,:)/dFnuOBS(x,y,:),MASK=mask(x,y,:))
  !       IF (medSovN(x,y) < 0._DP) medSovN(x,y) = 0._DP

  !     END IF
  !   END DO
  ! END DO

  !! Read the instrumental covariance matrix (TBD)

  !! Constraints
  ALLOCATE(itied(Npar))
  itied(:) = 0
  DO i=1,Npar
    IF (.NOT. TRIMEQ(parinfo(i)%tied,"")) &
      CALL IWHERE(TRIMEQ(parinfo(:)%name,parinfo(i)%tied),itied(i))
  END DO
  ! CALL IWHERE(( .NOT. parinfo(:)%fixed ) .AND. ( itied(:) == 0 ),ifree)
  ! Ncons = NwOBS - Nfreefilt
  ! Nfreepar = COUNT( (.NOT. parinfo(:)%fixed ) .AND. ( itied(:) == 0 ) )
  ! Ndof = Ncons - Nfreepar

  !! c. Compute the covariance matrix of the uncertainties
  !!-------------------------------------------------------
  !! Mask for the RMS covariance (diagonal)
  ALLOCATE(maskS(NwOBS,NwOBS))
  maskS(:,:) = .FALSE.
  FORALL (i=1:NwOBS) maskS(i,i) = .TRUE.

  !! 4D mask for the covariance matrix
  ! ALLOCATE(mask4D(Nx,Ny,NwOBS,NwOBS))
  ! mask4D(:,:,:,:) = .TRUE.
  ! FORALL (x=1:Nx,y=1:Ny,i=1:NwOBS, .NOT. mask(x,y,i))
  !   mask4D(x,y,i,:) = .FALSE.
  !   mask4D(x,y,:,i) = .FALSE.
  ! END FORALL

  !! Covariance matrix for each pixel and its inverse
  ALLOCATE(covarOBS(Nx,Ny,NwOBS,NwOBS), invLcovarOBS(Nx,Ny,NwOBS,NwOBS))

  covarOBS(:,:,:,:) = 0._DP
  !! invLcovarOBS is a lower triangle matrix L (Likelihood) results from
  !! the Cholesky decomposition of the covariance matrix (Appendix C2, Galliano18a)
  invLcovarOBS(:,:,:,:) = 0._DP
  FORALL (x=1:Nx,y=1:Ny,ALL(mask(x,y,:)))
    covarOBS(x,y,:,:) = UNPACK(dFnuOBS(x,y,:)**2, maskS, FIELD=0._DP)
    invLcovarOBS(x,y,:,:) = UNPACK(1._DP/dFnuOBS(x,y,:), maskS, FIELD=0._DP)

  END FORALL

  !! No wave calib (between diff instr/module)
  ! DO x=1,Nx
  !   DO y=1,Ny
  !     DO i=1,NwOBS
  !       IF (.NOT. mask(x,y,i)) THEN
  !         invLcovarOBS(x,y,i,:) = 0._DP
  !         invLcovarOBS(x,y,:,i) = 0._DP

  !       END IF
  !     END DO
  !   END DO
  ! END DO

  ! DEALLOCATE(covarOBS)

  print*, 'Read the inputs [done]'//NEW_LINE('')

  print*, '=================================================='

  !!------------------------------------------------------------------------
  !!                          II. Initial parameters
  !!------------------------------------------------------------------------

  DO i=1,MERGE(2,1,verbose)
    WRITE(unitlog(i),*)
    WRITE(unitlog(i),*) "ESTIMATING THE INITIAL VALUES OF THE PARAMETERS (" &
                        //TRIMLR(TIMINFO(timestr))//")"
  END DO

  !! a. Rough spectrum estimators
  !!------------------------------

  !! b. Automatic init param estimates
  !!-----------------------------------
  ALLOCATE(parini(Nx, Ny, Npar, MAX(NiniMC,1)))
  parini(:,:,:,:) = 0._DP
  CALL INITPARAM(NiniMC, IND=ind, PAR=parini(:,:,:,:), PARINFO=parinfo(:), &
                 ITIED=itied(:), MASK=mask(:,:,:), &
                 NEWINIT=newinit, FILOBS=filOBS, &
                 LABB=labB(:), LABL=labL(:), QABS=Qabs(:))

  Nparfree = COUNT((.NOT. parinfo(:)%fixed) .AND. (itied(:) <= 0))
  PRINT*, 'Number of free param: '//TRIMLR(PRING(Nparfree))

  print*, 'Initial parameters [done]'//NEW_LINE('')

  print*, '=================================================='

  !!------------------------------------------------------------------------
  !!                            Run the fitter
  !!------------------------------------------------------------------------

  ALLOCATE(par(Nx,Ny,Npar), status(Nx,Ny), resid(NwOBS), chi2red(Nx,Ny), &
           parerr(Nx,Ny,Npar), covpar(Nx,Ny,Npar,Npar), Niter(Nx,Ny), &
           Fnu_mod(NwOBS), limits(Npar,2), limited(Npar,2))

  FORALL (i=1:Npar)
    limits(i,:) = parinfo(i)%limits(:)
    limited(i,:) = parinfo(i)%limited(:)

  END FORALL

  ALLOCATE (statusMC(MAX(NiniMC,1)),NiterMC(MAX(NiniMC,1)), &
            parerrMC(Npar,MAX(NiniMC,1)),chi2redMC(MAX(NiniMC,1)), &
            covarMC(Npar,Npar,MAX(NiniMC,1)))

  par(:,:,:) = parini(:,:,:,1)
  CALL WRITE_HDF5(INITDBLARR=[Nx,Ny,Npar], NAME='Best fitted parameter value', &
                  FILE=filOUT, COMPRESS=compress, VERBOSE=debug, APPEND=.FALSE.)

  chi2fitx: DO xOBS=1,Nx
    chi2fity: DO yOBS=1,Ny
      notmasked: IF (ANY( mask(xOBS,yOBS,:) )) THEN
        PRINT*
        PRINT*, 'pos: ('//TRIMLR(PRING(xOBS))//', '//TRIMLR(PRING(yOBS))//'): '

        !! a. Levenberg-Marquardt method
        !!-------------------------------
        pariniMC: IF (NiniMC == 0) THEN

          CALL CHI2MIN_LM (residuals, NwOBS, PAR=par(xOBS,yOBS,:), VERBOSE=debug, &
                           STATUS=status(xOBS,yOBS), PARNAME=parinfo(:)%name, &
                           LIMITED=limited(:,:), LIMITS=limits(:,:), &
                           FIXED=parinfo(:)%fixed, ITIED=itied(:), &
                           CHI2RED=chi2red(xOBS,yOBS), NITER=Niter(xOBS,yOBS), &
                           PARERR=parerr(xOBS,yOBS,:), COVAR=covpar(xOBS,yOBS,:,:), NITERMAX=40)
        ELSE

          !! Vary the init param and keep the best fit
          DO j=1,NiniMC
            PRINT*, 'iniMC iteration: '//TRIMLR(PRING(j))//'/'//TRIMLR(PRING(NiniMC))
            CALL CHI2MIN_LM (residuals, NwOBS, PAR=parini(xOBS,yOBS,:,j), VERBOSE=debug, &
                             STATUS=statusMC(j), PARNAME=parinfo(:)%name, &
                             LIMITED=limited(:,:), LIMITS=limits(:,:), &
                             FIXED=parinfo(:)%fixed, ITIED=itied(:), &
                             CHI2RED=chi2redMC(j), NITER=NiterMC(j), &
                             PARERR=parerrMC(:,j), COVAR=covarMC(:,:,j), NITERMAX=40)

            DO i=1,MERGE(2,1,verbose)
              WRITE(unitlog(i),*)
              WRITE(unitlog(i),*) "PROGRAM EXECUTED IN "//TRIMLR(TIMINFO(timestr))//"."
              WRITE(unitlog(i),*)
            END DO

          END DO
          ibest = MAXVAL(MINLOC(chi2redMC(:)))
          par(xOBS,yOBS,:) = parini(xOBS,yOBS,:,ibest)
          parerr(xOBS,yOBS,:) = parerrMC(:,ibest)
          covpar(xOBS,yOBS,:,:) = covarMC(:,:,ibest)
          status(xOBS,yOBS) = statusMC(ibest)
          Niter(xOBS,yOBS) = NiterMC(ibest)
          chi2red(xOBS,yOBS) = chi2redMC(ibest)
          DO i=1,MERGE(2,1,verbose)
            WRITE(unitlog(i),*) "  MC chi2red is between " &
              //TRIMLR(PRING(MINVAL(chi2redMC(:)),NDEC=6))//" and " &
              //TRIMLR(PRING(MAXVAL(chi2redMC(:)),NDEC=6))

          END DO

        END IF pariniMC

        !! b. Derived quantities
        !!-----------------------

        !! c. Correlation between parameters
        !!-----------------------------------

        !! d. Summary
        !!------------
        IF (verbose) THEN
          PRINT*, '>>>'
          PRINT*, "Checking output ("//TRIMLR(PRING(xOBS))//', '//TRIMLR(PRING(yOBS))//'): '
          PRINT*, "  status = "//TRIMLR(PRING(status(xOBS,yOBS)))
          PRINT*, "  chi2red = "//TRIMLR(PRING(chi2red(xOBS,yOBS),NDEC=10))
          PRINT*, "  Niter = "//TRIMLR(PRING(Niter(xOBS,yOBS)))
          PRINT*, "Covariance matrix:"
          ! DO i=1,Npar
          !   PRINT*, REAL(covpar(xOBS,yOBS,:,i), KIND(0.))
          ! END DO
          PRINT*, '<<<'
          PRINT*

        END IF

        CALL WRITE_HDF5(DBLARR3D=par(xOBS:xOBS,yOBS:yOBS,:), &
                        NAME='Best fitted parameter value', &
                        FILE=filOUT, COMPRESS=compress, VERBOSE=debug, &
                        IND1=[xOBS,xOBS], IND2=[yOBS,yOBS])

      END IF notmasked
    END DO chi2fity
  END DO chi2fitx

  !! Extending init param for HB inputs
  CALL WRITE_HDF5(DBLARR3D=par(:,:,:), NAME='Chi2 fitted parameter value', &
                  FILE=filOBS, COMPRESS=compress, VERBOSE=debug, APPEND=.TRUE.)

  print*, 'Chi2 fit M83 IRS obs [done]'//NEW_LINE('')

  print*, '=================================================='

  !!------------------------------------------------------------------------
  !!                                 Analysis
  !!------------------------------------------------------------------------

  !! Calculate model
  !!-----------------
  ALLOCATE(FnuMOD(Nx,Ny,NwOBS))

  FnuMOD(:,:,:) = specModel( wOBS(:), INDPAR=ind, PARVAL=par(:,:,:), &
                             QABS=Qabs(:), EXTINCT=extinct(:,:), &
                             FNUCONT=FnuCONT, FNUBAND=FnuBAND, FNUSTAR=FnuSTAR, &
                             PABS=Pabs, FNULINE=FnuLINE, &
                             FNUCONT_TAB=FnuCONT_tab, FNUBAND_TAB=FnuBAND_tab, &
                             FNUSTAR_TAB=FnuSTAR_tab, PABS_TAB=Pabs_tab, &
                             FNULINE_TAB=FnuLINE_tab )

  DO i=1,Nextc
    FnuCONT_tab(:,:,:,i) = FnuCONT_tab(:,:,:,i) * Pabs(:,:,:)
  END DO
  CALL WRITE_HDF5(DBLARR4D=FnuCONT_tab, NAME='FnuCONT ('//TRIMLR(spec_unit)//')', &
                  FILE=filOUT, COMPRESS=compress, VERBOSE=debug, APPEND=.TRUE.)
  DO i=1,Nline
    FnuLINE_tab(:,:,:,i) = FnuLINE_tab(:,:,:,i) + (FnuCONT(:,:,:)+FnuSTAR(:,:,:))*Pabs(:,:,:)
  END DO
  CALL WRITE_HDF5(DBLARR4D=FnuLINE_tab, NAME='FnuLINE ('//TRIMLR(spec_unit)//')', &
                  FILE=filOUT, COMPRESS=compress, VERBOSE=debug, APPEND=.TRUE.)
  DO i=1,Nband
    FnuBAND_tab(:,:,:,i) = FnuBAND_tab(:,:,:,i) + (FnuCONT(:,:,:)+FnuSTAR(:,:,:))*Pabs(:,:,:)
  END DO
  CALL WRITE_HDF5(DBLARR4D=FnuBAND_tab, NAME='FnuBAND ('//TRIMLR(spec_unit)//')', &
                  FILE=filOUT, COMPRESS=compress, VERBOSE=debug, APPEND=.TRUE.)
  DO i=1,Nstar
    FnuSTAR_tab(:,:,:,i) = FnuSTAR_tab(:,:,:,i) * Pabs(:,:,:)
  END DO
  CALL WRITE_HDF5(DBLARR4D=FnuSTAR_tab, NAME='FnuSTAR ('//TRIMLR(spec_unit)//')', &
                  FILE=filOUT, COMPRESS=compress, VERBOSE=debug, APPEND=.TRUE.)
  CALL WRITE_HDF5(DBLARR3D=FnuMOD, NAME='FnuMOD ('//TRIMLR(spec_unit)//')', &
                  FILE=filOUT, COMPRESS=compress, VERBOSE=debug, APPEND=.TRUE.)

  print*, 'fitChi2 analysis [done]'//NEW_LINE('')

  print*, '=================================================='

  !! Final
  !!-------
  DO i=1,MERGE(2,1,verbose)
    WRITE(unitlog(i),*)
    WRITE(unitlog(i),*) "PROGRAM EXECUTED IN "//TRIMLR(TIMINFO(timestr))//"."
    WRITE(unitlog(i),*)
  END DO
  IF (verbose) THEN
    PRINT*, " - File "//fiLOG//" has been written."
    PRINT*
  END IF

  !! Free memory space
  DEALLOCATE(wOBS, FnuOBS, dFnuOBS, maskint, mask)

END PROGRAM test_fitChi2

max_stars_repo_name: kxxdhdn/MISSILE

max_issues_repo_name: kxxdhdn/MISSILE

max_forks_repo_name: kxxdhdn/MISSILE


!=======================================================================
! Occam's Inversion Fortran2003 Code
!=======================================================================
!
!    Copyright (C) 1986-2011
!    Steven Constable, Kerry Key, David Myer, Catherine deGroot-Hedlin
!    Scripps Institution of Oceanography
!    University of California, San Diego
!
!    This file is part of Occam's Inversion.
!
!    Occam's Inversion is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    Occam's Inversion is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with Occam's Inversion.  If not, see <http://www.gnu.org/licenses/>.
!
!-----------------------------------------------------------------------
!
! This file is a routine to run Occam's inversion.
! It calls module Occam in file Occam.f90 for computing
! each Occam iteration and for performing Occam I/O.
!
! For extensive Occam programming notes, see top of file Occam.f90.
!
!-----------------------------------------------------------------------
    program runOccam
!-----------------------------------------------------------------------

    ! Modules to include:
    use Occam

    implicit none

! local variables:
    integer  maxitr, itmax, konv, ifftol
!   lerr = error flag for file operations
!   maxitr = maximum number of subsequent iterations by this run
!   itmax = maximum iteration number
!   konv = convergence flag from OCCAM
!   ifftol = convergence info from OCCAM (has the desired misfit been achieved?)
    real(RealPrec) tolreq, tobt, pmu, stepsz, rlast
!   tolreq = desired misfit, in RMS standard errors
!   tobt = current misfit
!   pmu = current Lagrange multiplier
!   stepsz = current RMS step through model space
!   rlast = last RMS measure of roughness
    character(50) descr,  cStartup
!   descr = description line in iteration files
!   cStartup = name of starting file.  (See notes below)

    logical   :: bFwdOnly
!   bFwdOnly = set true if -F on cmd line ... saying to produce a fwd model & exit.
    logical :: lrunOccam
!   lrunOccam = .true. to run the occam iterations.  Some compilers do not
!     like looping over a static .true.

    real(8) :: timeOffset ! Timer temp:
    real(8) :: time0      ! Start time of Occam program (in seconds)
    real(8) :: timeIter0  ! Start time of an Occam iteration

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

!
! Start the timer:
!
    call get_time_offset(0d0, time0)

! DGM Oct 2006 - look for a command line parameter telling the name of
! the startup file.  If not given, then use 'startup'.
! If two params are given, then first is the startup file and the 2nd
! is the root filename for output of iteration, response, & log files.
    write(*,*) ' '
    write(*,*) '------------------------------------------------------------'
    write(*,*) ' Occam''s Inversion                Version 3.10, March 22, 2011'
    write(*,*) ' '
    write(*,*) ' Original Reference: '
    write(*,*) ''
    write(*,*) ' Constable, S. C., R. L. Parker, and C. G. Constable, 1987, '
    write(*,*) ' Occam''s inversion - A practical algorithm for generating   '
    write(*,*) ' smooth models from electromagnetic sounding data, '
    write(*,*) ' Geophysics, 52, 289-300. '
    write(*,*) ' '
    write(*,*) ' Copyright (C) 1986-2011'
    write(*,*) ' Steven Constable, Kerry Key, David Myer, Catherine deGroot-Hedlin'
    write(*,*) ' Scripps Institution of Oceanography'
    write(*,*) ' University of California, San Diego'
    write(*,*) ' '
    write(*,*) ' License: '
    write(*,*) ' '
    write(*,*) ' Occam''s Inversion is free software: you can redistribute it and/or modify'
    write(*,*) ' it under the terms of the GNU General Public License as published by'
    write(*,*) ' the Free Software Foundation, either version 3 of the License, or'
    write(*,*) ' (at your option) any later version.'
    write(*,*) ' '
    write(*,*) ' Occam''s Inversion is distributed in the hope that it will be useful,'
    write(*,*) ' but WITHOUT ANY WARRANTY; without even the implied warranty of'
    write(*,*) ' MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the'
    write(*,*) ' GNU General Public License for more details.'
    write(*,*) '  '
    write(*,*) ' You should have received a copy of the GNU General Public License'
    write(*,*) ' along with Occam''s Inversion.  If not, see <http://www.gnu.org/licenses/>.'
    write(*,*) ' '
    write(*,*) '------------------------------------------------------------'
    write(*,*) ' '
    cStartup    = 'startup'     ! Establish defaults
    cRootName   = 'ITER'        ! Default output file name root:
    bFwdOnly    = .false.

    !
    ! KWK, Dec 2008:
    ! New code for getting command line arguments using Fortran2003 standard intrinsics
    ! command_argument_count()  and get_command_argument()
    !
    call getOccamCommandLineArguments(cStartup,cRootName,bFwdOnly)

    ! Read in the startup file:
    ! (Startup and iteration files should be machine written, and so we do
    ! not necessarily test the validity of the input data unless arrays are
    ! involved.)
    call readIteration (cStartup,descr,maxitr,tolreq, pmu,rlast,tobt,ifftol)

    ! Read in the data file:
    call readData

    ! Allocate things related to the number of data
    call allocateNdArrays

    ! Read in the model:
    call readModel

    ! Allocate working arrays for storing associated parameters if defined:
    call allocateAssocArrays

    ! If we are to produce a forward model only, do so & exit.
    if (bFwdOnly) then
        write(*,*) 'Outputting FORWARD MODEL ONLY...'

        call computeFwd( .false., pm )

        ! Show the misfit
        rlast = sqrt(sum( ((D - DM)/SD)**2 )/nd )       ! Array math - stack alloc of only nd size

        write(*,*) ' Forward Model Response Misfit:', rlast

        ! write out the model response:
        if (trim(cRootName) == 'ITER') then
            ! Assign a default rootname:
            cRootName = cStartup ! just point to the startupfile
        endif
        call writeResponse()
        write(*,*) ' '
        write(*,*) ' Done! '
        write(*,*) ' '

    else        ! Invert!

        ! Open a logfile:
        call openOccamLogFile()

        ! if this is zeroth iteration copy the startup file into ITER00
        if (nCurrentIter == 0) then
            ! initialize the inversion parameters
            pmu     = 5.0
            rlast   = 1.0e+10
            tobt    = 1000.0
            ifftol  = 0
            call writeIteration (descr,maxitr,tolreq,pmu,rlast,tobt,ifftol)
        end if

        itmax = maxitr + nCurrentIter
        write(*,*) 'Major I/O accomplished, sit back and relax.......'
        write(*,*) ' '

        lrunOccam= .true.

        ! Run the main loop
        do while (lrunOccam)

            ! Start the iteration timer:
            call get_time_offset( 0d0,timeIter0)

            ! Compute an Occam iteration:
            call compOccamIteration(tolreq,tobt,ifftol,pmu,rlast,stepsz,konv)

            ! DGM Feb 2010 if the convergence code is 2 or 3, then the inversion
            ! is ending without converging.  In this case, dm (model response)
            ! is not guaranteed to be a realistic model.  It is sometimes set
            ! to something wild from the various minimizer routines searching
            ! desperately for a solution. Do NOT output this iteration & response
            ! because it is worthless anyway.  In these cases the rms misfit reported
            ! in the iteration file does not match what you can calculate from
            ! the response file.
            if (.not.(konv == 2 .or. konv == 3)) then
                ! write an iteration output file:
                call writeIteration ( descr,maxitr,tolreq,pmu,rlast,tobt,ifftol )

                ! write out the model response:
                call writeResponse()
            endif

            ! Display information to the terminal window:
            write(*,*) ' '
            write(*,*) ' ----------------------------------'
            write(*, '(a24,i8)') ' Iteration:', nCurrentIter
            write(*, '(a24,g12.4)') ' RMS Misfit:', tobt
            if (ifftol == 1) then
                write(*, '(a24,g12.4,a9)') 'Target Misfit:', tolreq, ' reached!'
            else
                write(*, '(a24,g12.4)') 'Target Misfit:', tolreq
            endif
            write(*, '(a24,g12.4)') ' Roughness:', rlast
            write(*, '(a24,g12.4)') ' Stepsize:', stepsz
            call get_time_offset( timeIter0,timeOffset)
            write(*, '(a24,g12.4)') ' Iteration Time (s):', timeOffset
            call get_time_offset( time0,timeOffset)
            write(*, '(a24,g12.4)') ' Total Time (s):', timeOffset
            write(*,*) ' ----------------------------------'
            write(*,*) ' '

            ! test for convergence: continue inversion if
            !   1) iteration number does not exceed maximum, and
            !   2) stepsize is still significant, and
            !   3) rms misfit larger than required, and
            !   4) there are no irrecoverable convergence problems
            if (konv .ne. 0) then
                select case (konv)
                case (1)
                    write(ioUnitOccamLogFile,*)   'Perfectly smooth model found.'
                    write(*,*)      'Perfectly smooth model found.'
                case (2)
                    write(ioUnitOccamLogFile,*)   'Convergence problems from RMS Misfit'
                    write(*,*)      'Convergence problems from RMS Misfit'
                    write(*,*)      'May have reached minimum possible RMS Misfit level!'
                case (3)
                    write(ioUnitOccamLogFile,*)   'Convergence problems from Smoothness'
                    write(*,*)      'Convergence problems from Smoothness'
                    write(*,*)      'May have reached minimum possible RMS Misfit level!'
                case (4)
                    write(ioUnitOccamLogFile,*)   'Misfit changed by less than Delta Misfit Limit.'
                    write(*,*)      'Misfit changed by less than Delta Misfit Limit.'
                    write(*,*)      'Stopping inversion so you can check the data to see'
                    write(*,*)      'if the error level needs to be adjusted or if some'
                    write(*,*)      'data points need to be removed.'
                case default
                    write(ioUnitOccamLogFile,*)   ' Convergence problems of type ', konv
                    write(*,*)      ' Convergence problems of type ', konv
                end select
                close(ioUnitOccamLogFile)
                exit

            else if (nCurrentIter >= itmax) then
            ! Inversion has reached max # of iterations specified in startup file.
            !   Stop now even if not converged.
            ! NOTE: This check MUST be before the "continue" check below.
                write(ioUnitOccamLogFile,*) ' Stop on maximum iterations'
                close(ioUnitOccamLogFile)
                write(*,*) ' Reached maximum iterations.  Done.'
                exit

            else if ((tobt < 1.01*tolreq) .and. (stepsz < 0.001)) then
                write(ioUnitOccamLogFile,*) ' Stop on normal convergence'
                close(ioUnitOccamLogFile)
                write(*,*) ' Solution converged!  Done.'
                exit
            !
            ! else... continue...

            end if
        enddo
    endif

    ! Deallocate memory & return control to the user
    call deallocateOccam

end program runOccam


text:
+incdir+../testbench
+incdir+../cores/orpsoc-cores/cores/wb_bfm/wb_bfm-0/
+incdir+../cores/orpsoc-cores/cores/verilog_utils/
+define+VERBOSE
+define+SIM

../rtl/syscon.v
../rtl/testing_wb_slave.v
../testbench/testbench_testing_wb_slave.v
../testbench/arbiter.v
../testbench/wb_switch1.v
../testbench/dump.v
../cores/orpsoc-cores/cores/wb_bfm/wb_bfm-0/wb_bfm_master.v
../cores/orpsoc-cores/cores/wb_intercon/wb_mux.v
../cores/orpsoc-cores/cores/wb_intercon/wb_arbiter.v
../cores/orpsoc-cores/cores/wb_intercon/wb_data_resize.v


text:
A bike circle is an easy way for many Bike bikes to cross without anyone killing each other. At least, thats how its supposed to work. The Wiki Logo Davis Wiki Logo is based on the bike circle signs.

Always go around them counterclockwise (as seen from above). If you dont, someones bound to spot you and think youre a freshman. Worse yet, you could cause an accident. If youre a pedestrian, USE THE CROSSWALKS and dont walk through the bike circle. It slows traffic down and increases the risk of you getting hit.

Also critical to smooth bikecircling is to observe the proper rightofway: Bikes already in the circle have the rightofway over bikes entering the circle.

The area in the middle of a bike circle is usually thought of a Nomansland, since no one bikes through it, and few people rarely walk across it. It would be a perfect spot for a picnic in the fall, since youd have a good view of all the bike crashes  bring lawn chairs.

Similar schemes called Traffic Circle traffic circles or roundabouts, meant to accommodate motor traffic, are in place on roads in Davis and West Village.

In late August 2007, some of the Aggie Roundabout bike circles were painted with Aggie logos.

Pictures

20051011 10:59:23 nbsp Putting those large concrete bumps there, and on bike paths, is hazardous the cure is far worse than the problem. They should have just painted lines, not the bricks. I think there should be a movement to remove the bumps. Users/SteveDavison

20051011 11:03:15 nbsp Since the space inside is wasted, why not fill it with a large round flower planter? It would be more even more effective than the bumps, plus it would add beauty and be environmentally friendly. Intersections waste a lot of space, enough for a large house look at the picture. Users/SteveDavison
    Part of the reason the center area isnt filled is to allow large vehicles to use the intersection
    Because me and my friends get a kick out of setting up a few lawnchairs in the center with a few cold drinks and watching the newbies scream unneccessarily, crash, and start a domino effect around the entire circle.
    You also dont want anything in the middle of the circle so that you can easily see other bikes approaching the circle. Users/JimEvans

20060525 18:51:24 nbsp I actually for the most part WILL ride my bike right on through those bike circles(going in circles at 20mph is a pain(I wiped out once that way)) Users/DevinCastellucci

20061029 12:07:02 nbsp Part of the point of the circles at the intersection is to make you slow down. In addition to slower speeds being safer, it also makes all the traffic flow in one direction, which is safer than your standard intersection.  Users/MatthewTom

It really sucks when a large pack of roving students enters in the circle at the same time, forcing you to flow with them or cause a major incident, people go to slowly anyway. ~Users/Davidpoole

20071015 23:04:44 nbsp Yield and use turn signals. For the love of god, yield and use turn signals! Users/RobertSooHoo

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search