get_fio_GCmagnetic_fields_p Subroutine

  subroutine get_fio_GCmagnetic_fields_p(params,F,Y_R,Y_PHI,Y_Z, &
       B_R,B_PHI,B_Z,gradB_R,gradB_PHI,gradB_Z, &
       curlb_R,curlb_PHI,curlb_Z,flag,hint)
    TYPE(FIELDS), INTENT(IN)       :: F
    TYPE(KORC_PARAMS), INTENT(IN)  :: params
    REAL(rp), DIMENSION(params%pchunk), INTENT(IN)  :: Y_R,Y_PHI,Y_Z
    REAL(rp), DIMENSION(params%pchunk), INTENT(INOUT)  :: B_R,B_PHI,B_Z
    REAL(rp), DIMENSION(params%pchunk), INTENT(INOUT)  :: gradB_R,gradB_PHI,gradB_Z
    REAL(rp), DIMENSION(params%pchunk), INTENT(INOUT)  :: curlb_R,curlb_PHI,curlb_Z
    INTEGER(is), DIMENSION(params%pchunk), INTENT(INOUT)  :: flag
    TYPE(C_PTR), DIMENSION(params%pchunk), INTENT(INOUT)  :: hint
    INTEGER (C_INT)                :: status
    INTEGER                        :: pp,pchunk
    REAL(rp), DIMENSION(3)         :: x
    REAL(rp), DIMENSION(3)         :: Btmp
    REAL(rp), DIMENSION(9)         :: dBtmp
    REAL(rp)  :: Bmag,dBRdR,dBPHIdR,dBZdR,dBRdPHI,dBPHIdPHI,dBZdPHI,dBRdZ,dBPHIdZ,dBZdZ

    pchunk=params%pchunk

    do pp = 1,pchunk
       if (flag(pp) .EQ. 1_is) then
          x(1) = Y_R(pp)*params%cpp%length
          x(2) = Y_PHI(pp)
          x(3) = Y_Z(pp)*params%cpp%length

          !             prtcls%hint(pp)=c_null_ptr

          !if (pp.eq.1) write(output_unit_write,*) 'Yinterp',x

          status = fio_eval_field(F%FIO_B, x(1), &
               Btmp(1),hint(pp))

          if (status .eq. FIO_SUCCESS) then

             !if (pp.eq.1) write(output_unit_write,*) 'interp success!'
             
             B_R(pp)=Btmp(1)/params%cpp%Bo
             B_PHI(pp)=Btmp(2)/params%cpp%Bo
             B_Z(pp)=Btmp(3)/params%cpp%Bo

             Bmag=sqrt(B_R(pp)*B_R(pp)+B_PHI(pp)*B_PHI(pp)+B_Z(pp)*B_Z(pp))
             
             status = fio_eval_field_deriv(F%FIO_B, x(1),dBtmp(1),hint(pp))

             !dBRdR=dBtmp(FIO_DR_R)*(params%cpp%length/params%cpp%Bo)
             !dBPHIdR=dBtmp(FIO_DR_PHI)*(params%cpp%length/params%cpp%Bo)
             !dBZdR=dBtmp(FIO_DR_Z)*(params%cpp%length/params%cpp%Bo)
             !dBRdPHI=dBtmp(FIO_DPHI_R)*(params%cpp%length/params%cpp%Bo)
             !dBPHIdPHI=dBtmp(FIO_DPHI_PHI)*(params%cpp%length/params%cpp%Bo)
             !dBZdPHI=dBtmp(FIO_DPHI_Z)*(params%cpp%length/params%cpp%Bo)
             !dBRdZ=dBtmp(FIO_DZ_R)*(params%cpp%length/params%cpp%Bo)
             !dBPHIdZ=dBtmp(FIO_DZ_PHI)*(params%cpp%length/params%cpp%Bo)
             !dBZdZ=dBtmp(FIO_DZ_Z)*(params%cpp%length/params%cpp%Bo)

             dBRdR=dBtmp(1)*(params%cpp%length/params%cpp%Bo)
             dBPHIdR=dBtmp(2)*(params%cpp%length/params%cpp%Bo)
             dBZdR=dBtmp(3)*(params%cpp%length/params%cpp%Bo)
             dBRdPHI=dBtmp(4)/params%cpp%Bo
             dBPHIdPHI=dBtmp(5)/params%cpp%Bo
             dBZdPHI=dBtmp(6)/params%cpp%Bo
             dBRdZ=dBtmp(7)*(params%cpp%length/params%cpp%Bo)
             dBPHIdZ=dBtmp(8)*(params%cpp%length/params%cpp%Bo)
             dBZdZ=dBtmp(9)*(params%cpp%length/params%cpp%Bo)

             !write(6,*) 'dBRdR',dBRdR
             !write(6,*) 'dBPHIdR',dBPHIdR
             !write(6,*) 'dBZdR',dBZdR
             !write(6,*) 'dBRdPHI',dBRdPHI
             !write(6,*) 'dBPHIdPHI',dBPHIdPHI
             !write(6,*) 'dBZdPHI',dBZdPHI
             !write(6,*) 'dBRdZ',dBRdZ
             !write(6,*) 'dBPHIdZ',dBPHIdZ
             !write(6,*) 'dBZdZ',dBZdZ
             
             gradB_R(pp)=(B_R(pp)*dBRdR+B_PHI(pp)*dBPHIdR+B_Z(pp)*dBZdR)/ &
                  Bmag
             gradB_PHI(pp)=(B_R(pp)*dBRdPHI+B_PHI(pp)*dBPHIdPHI+ &
                  B_Z(pp)*dBZdPHI)/(Y_R(pp)*Bmag)
             gradB_Z(pp)=(B_R(pp)*dBRdZ+B_PHI(pp)*dBPHIdZ+B_Z(pp)*dBZdZ)/ &
                  Bmag

             curlb_R(pp)=(Bmag/Y_R(pp)*dBZdPHI-B_Z(pp)*gradB_PHI(pp)- &
                  Bmag*dBPHIdZ+B_PHI(pp)*gradB_Z(pp))/(Bmag*Bmag)
             curlb_PHI(pp)=(Bmag*dBRdZ-B_R(pp)*gradB_Z(pp)- &
                  Bmag*dBZdR+B_Z(pp)*gradB_R(pp))/(Bmag*Bmag)
             curlb_Z(pp)=(Bmag/Y_R(pp)*B_PHI(pp)+Bmag*dBPHIdR- &
                  B_PHI(pp)*gradB_R(pp)- &
                  Bmag/Y_R(pp)*dBRdPHI+B_R(pp)*gradB_PHI(pp))/(Bmag*Bmag)
             
          else if (status .eq. FIO_NO_DATA) then            
             flag(pp) = 0_is
          else if (status .ne. FIO_SUCCESS) then
             flag(pp) = 0_is
          end if
          
       end if
    end do

  end subroutine get_fio_GCmagnetic_fields_p