diff --git a/debian/changelog b/debian/changelog
index 7ae0ae3..986adea 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,9 @@
+healpix-fortran (3.81-1) UNRELEASED; urgency=low
+
+  * New upstream release.
+
+ -- Debian Janitor <janitor@jelmer.uk>  Wed, 27 Apr 2022 00:40:46 -0000
+
 healpix-fortran (3.60+ds-1) unstable; urgency=medium
 
   * Initial release. (Closes: #963092)
diff --git a/debian/patches/0001-Link-against-libsharp-and-libcfitsio.patch b/debian/patches/0001-Link-against-libsharp-and-libcfitsio.patch
index a9ba311..baab267 100644
--- a/debian/patches/0001-Link-against-libsharp-and-libcfitsio.patch
+++ b/debian/patches/0001-Link-against-libsharp-and-libcfitsio.patch
@@ -6,10 +6,10 @@ Subject: Link against libsharp and libcfitsio
  src/f90/mod/Makefile | 2 +-
  1 file changed, 1 insertion(+), 1 deletion(-)
 
-diff --git a/src/f90/mod/Makefile b/src/f90/mod/Makefile
-index d0a45f3..2344ed7 100644
---- a/src/f90/mod/Makefile
-+++ b/src/f90/mod/Makefile
+Index: healpix-fortran/src/f90/mod/Makefile
+===================================================================
+--- healpix-fortran.orig/src/f90/mod/Makefile
++++ healpix-fortran/src/f90/mod/Makefile
 @@ -51,7 +51,7 @@ endif
  all: $(OD)_mkdir $(LIB)
  
diff --git a/src/._healpy b/src/._healpy
new file mode 100755
index 0000000..4dd22d6
Binary files /dev/null and b/src/._healpy differ
diff --git a/src/f90/lib/gifmod.f90 b/src/f90/lib/gifmod.f90
index a30276f..a9feeda 100644
--- a/src/f90/lib/gifmod.f90
+++ b/src/f90/lib/gifmod.f90
@@ -222,7 +222,7 @@ module gifmod
      subroutine gifstr(x,nx,ny,sx,sy,nc,r,g,b,or_,fn) bind(C, name="gifstr")
      use iso_c_binding
      implicit none
-     integer (c_int), intent(out) :: x(*)
+     integer (c_int), intent(inout) :: x(*) ! 2020-08-25
      integer (c_int), intent(in) :: nx, ny, sx, sy, nc, r(*), g(*), b(*), or_
      character(c_char), intent(in) :: fn(*)
      end subroutine
@@ -485,8 +485,9 @@ contains
    end subroutine addttl
 
    ! ------------------------------------------------------------------
-   ! addstr: used exclusively by addbar to add labels to the ends of
-   !         color bar; [a] is the image array, [i] and [j] are the
+   ! addstr: used by addbar to add labels to the ends of color bar 
+   !          and by addttl to add title; 
+   !         [a] is the (modified) image array, [i] and [j] are the
    !         array positions where the string [s] should end if [or]
    !         is negative, or start if [or] is positive
    ! ------------------------------------------------------------------
@@ -494,10 +495,11 @@ contains
    subroutine addstr(a,i,j,or,s)
       use healpix_types
       implicit none
-      integer(i4b) :: nx,ny,i,j
+      integer(i4b), intent(inout), dimension(:,:) :: a
+      integer(i4b), intent(in) :: i,j
       integer(i4b), intent(in) :: or
-      integer(i4b), dimension(:,:) :: a
       character(len=*), intent(in) :: s
+      integer(i4b) :: nx,ny
 
       nx = size(a,1)
       ny = size(a,2)
diff --git a/src/f90/mod/Makefile b/src/f90/mod/Makefile
index d0a45f3..b3985e3 100644
--- a/src/f90/mod/Makefile
+++ b/src/f90/mod/Makefile
@@ -13,7 +13,7 @@ incpix = convert_inplace_1d_inc.f90 \
 	apply_mask_inc.f90 \
 	pixel_routines.F90
 
-incfits = fits_s_inc.f90 fits_d_inc.f90
+incfits = fits_s_inc.F90 fits_d_inc.F90
 
 incudgrade = udgrade_s_inc.f90 udgrade_d_inc.f90
 
@@ -109,7 +109,7 @@ $(OD)/alm_tools.o:alm_tools.F90  $(incalm)
 
 $(OD)/fitstools.o: fitstools.F90 $(incfits)
 
-$(OD)/udgrade_nr.o: udgrade_nr.f90 $(incudgrade)
+$(OD)/udgrade_nr.o: udgrade_nr.F90 $(incudgrade)
 
 $(OD)/num_rec.o: num_rec.F90 $(incnumrec)
 
@@ -123,10 +123,10 @@ alm_map_ss_inc.F90: alm_map_template.F90 gen_alm_code
 alm_map_dd_inc.F90: alm_map_template.F90 gen_alm_code alm_map_ss_inc.F90
 	@./gen_alm_code
 
-fits_s_inc.f90: fits_template.f90 gen_fits_code
+fits_s_inc.F90: fits_template.F90 gen_fits_code
 	@./gen_fits_code
 
-fits_d_inc.f90: fits_template.f90 gen_fits_code fits_s_inc.f90
+fits_d_inc.F90: fits_template.F90 gen_fits_code fits_s_inc.F90
 	@./gen_fits_code
 
 udgrade_s_inc.f90: udgrade_template.f90 gen_udgrade_code
diff --git a/src/f90/mod/alm_map_template.F90 b/src/f90/mod/alm_map_template.F90
index 3fc2921..239d76b 100644
--- a/src/f90/mod/alm_map_template.F90
+++ b/src/f90/mod/alm_map_template.F90
@@ -327,7 +327,7 @@
 
     real(DP), dimension(1:2)         :: zbounds_in
     integer(I4B) :: l, m, ith                          ! alm related
-    integer(I8B) :: istart_south, istart_north, npix   ! map related
+    integer(I8B) :: istart_south, istart_north, npix, p! map related
     integer(I4B) :: nrings, nphmx
 
     real(DP),     dimension(1:4)              :: b_even, b_odd
@@ -354,13 +354,27 @@
 #ifdef USE_SHARP
     zbounds_in = (/-1.d0 , 1.d0/)
     if (present(zbounds)) zbounds_in = zbounds
-
+!     aspin = abs(spin)
+!     if ((aspin>0).and.(aspin<=100)) then
+!       call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
+!         alm(1:2,0:nlmax,0:nmmax),map(0:12*nsmax*nsmax-1,1:2), zbounds_in)
+!       return
+!     endif
+    ! 2021-06-15: use scalar routine when s=0
     aspin = abs(spin)
-    if ((aspin>0).and.(aspin<=100)) then
-      call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
-        alm(1:2,0:nlmax,0:nmmax),map(0:12*nsmax*nsmax-1,1:2), zbounds_in)
-      return
+    npix  = (12_I8B*nsmax)*nsmax
+    if (aspin > 0) then
+       call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
+            alm(1:2,0:nlmax,0:nmmax),map(0:npix-1,1:2), zbounds_in)
+    else
+       call sharp_hp_alm2map_x_KLOAD(nsmax,nlmax,nmmax,&
+            alm(1:1,0:nlmax,0:nmmax),map(0:npix-1,1), zbounds_in)
+       do p = 0_I8B, npix - 1_I8B
+          map(p,1) = -map(p,1) ! mimic alm sign convention of alm2map_spin
+          map(p,2) = 0.0_KMAP  ! unused for s=0
+       enddo
     endif
+    return
 #endif
     !=======================================================================
 
@@ -2420,7 +2434,6 @@
     if (present(zbounds)) zbounds_in = zbounds
     w8ring_in  = 1.d0
     if (present(w8ring))  w8ring_in  = w8ring
-
     call sharp_hp_map2alm_x_KLOAD(nsmax,nlmax,nmmax,map,alm,zbounds_in,w8ring_in)
 #else
 
@@ -2696,12 +2709,29 @@
     if (present(w8ring))  w8ring_in  = w8ring
 
 #ifdef USE_SHARP
+!     aspin = abs(spin)
+!     if ((aspin>0).and.(aspin<=100)) then
+!       call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
+!         map(0:12*nsmax*nsmax-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
+!       return
+!     endif
+    ! 2021-06-15: use scalar routine when s=0
     aspin = abs(spin)
-    if ((aspin>0).and.(aspin<=100)) then
-      call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
-        map(0:12*nsmax*nsmax-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
-      return
+    npix  = (12_I8B*nsmax)*nsmax
+    if (aspin > 0) then
+       call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
+            map(0:npix-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
+    else
+       call sharp_hp_map2alm_x_KLOAD(nsmax,nlmax,nmmax, &
+            map(0:npix-1,1),alm(1:1,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
+       do l=0,nlmax
+          do m=l,nmmax
+             alm(1,l,m) = -alm(1,l,m) ! implement alm sign convention of map2alm_spin
+             alm(2,l,m) = 0.0_DPC     ! unused for s=0
+          enddo
+       enddo
     endif
+    return
 #endif
 
     ! Healpix definitions
diff --git a/src/f90/mod/fits_template.f90 b/src/f90/mod/fits_template.F90
similarity index 76%
rename from src/f90/mod/fits_template.f90
rename to src/f90/mod/fits_template.F90
index 337094a..f5dab00 100644
--- a/src/f90/mod/fits_template.f90
+++ b/src/f90/mod/fits_template.F90
@@ -39,6 +39,8 @@
 !    subroutine write_bintabh_KLOAD         (8)
 !    subroutine unfold_weights_KLOAD
 !    subroutine unfold_weightslist_KLOAD
+!    subroutine read_fits_partial_KLOAD     (4/8)
+!    subroutine write_fits_partial_KLOAD     (4)
 !
 !
 ! K M A P   : map kind                 either SP or DP
@@ -175,19 +177,23 @@ subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
 
     LOGICAL(LGT) :: anynull, do_polcconv
     ! Note : read_fits_cut4 still SP and I4B only
-    integer(I4B), dimension(:), allocatable :: pixel
-!     real(KMAP),     dimension(:), allocatable :: signal
-    real(SP),     dimension(:), allocatable :: signal
+    integer(I4B), dimension(:),   allocatable :: pixel4
+    real(SP),     dimension(:),   allocatable :: signal
+    ! while read_fits_partial is SP/DP and I4B/I8B
+    integer(I8B), dimension(:),   allocatable :: pixel8
+    real(KMAP),   dimension(:,:), allocatable :: iqu
     integer(I4B) :: status
     integer(I4B) :: type_fits, nmaps_fits
     CHARACTER(LEN=80)  :: units1
     CHARACTER(LEN=80),dimension(1:30)  :: unitsm
     character(len=80), dimension(:), allocatable  :: hbuffer
 !    CHARACTER(LEN=80),dimension(:), allocatable  :: unitsm
-    integer(I4B) :: extno_i, extno_f
+    integer(I4B) :: extno_i, extno_f, n_ext
     CHARACTER(LEN=80)  :: polcconv
     integer(I4B)       :: ipolcconv, polar, nside_fits
     character(len=*), parameter :: primer_url = 'http://healpix.sf.net/pdf/intro.pdf'
+    character(len=*), parameter :: code = 'input_map'
+    
     !-----------------------------------------------------------------------
 
     units1    = ' '
@@ -208,6 +214,7 @@ subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
        do_polcconv = .not. ignore_polcconv
     endif
 
+    n_ext = getnumext_fits(filename)
     obs_npix = getsize_fits(filename, nmaps = nmaps_fits, type=type_fits, extno=extno_i, &
          polcconv=ipolcconv, polarisation=polar, nside=nside_fits)
 
@@ -228,47 +235,70 @@ subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
           units(1:size(units)) = unitsm(1:size(units))
        endif
        call map_bad_pixels(map, fmissing, fmiss_effct, imissing, verbose=.false.)
-!        do imap = 1, nmaps
-!           anynull = .true.
-!           if (anynull) then
-!              imissing = 0
-!              do i=0,npixtot-1
-!                 if ( ABS(map(i,imap)/fmissing -1.) < 1.e-5 ) then
-!                    map(i,imap) = fmiss_effct
-!                    imissing = imissing + 1
-!                 endif
-!              enddo
-!           endif
-!        enddo
        if (imissing > 0) write(*,'(a,1pe11.4)') 'blank value : ' ,fmissing
-    else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3)) then
-       ! cut sky FITS file, reading 1 map (I) or 3 maps (I,Q,U), each from a different extension
+    else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3) &
+         & .and. (nmaps_fits > nmaps) .and. n_ext == 1) then
+       ! partial FITS file, reading 1 map (I) or 3 maps (I,Q,U) from the same extension
+       obs_npix = getsize_fits(filename)
+       allocate(pixel8(0:obs_npix-1),          stat = status)
+       call assert_alloc(status, code, "pixel8")
+       allocate(iqu   (0:obs_npix-1, 1:nmaps), stat = status)
+       call assert_alloc(status, code, "iqu")
+       call read_fits_partial(filename, pixel8, iqu, header=hbuffer)
+       if (present(header)) then
+          do i=1,nlheader
+             header(i) = hbuffer(i)
+          enddo
+       endif
+       if (present(units)) then
+          do imap=1,min(size(units),nmaps-1)
+             call get_card(hbuffer,'TUNIT'//string(imap+1,format='(i0.0)'), units(imap))
+          enddo
+       endif
+       maxpix = maxval(pixel8)
+       minpix = maxval(pixel8)
+       if (maxpix > (npixtot-1) .or. minpix < 0) then
+          print*,'map constructed from file '//trim(filename)//', with pixels in ',minpix,maxpix
+          print*,' wont fit in array with ',npixtot,' elements'
+          call fatal_error
+       endif
+       map(:,:)        = fmiss_effct
+       do imap=1,nmaps
+          map(pixel8(:),imap) = iqu(:,imap)
+       enddo
+       imissing = npixtot - obs_npix
+       deallocate(iqu)
+       deallocate(pixel8)
+!     else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3) .and. n_ext == nmaps) then
+!        ! cut sky FITS file, reading 1 map (I) or 3 maps (I,Q,U), each from a different extension
+    else if (type_fits == 3 .and. (nmaps == 3) .and. n_ext == nmaps) then
+       ! cut sky FITS file, reading 3 maps (I,Q,U), each from a different extension
        do imap = 1, nmaps
           extno_f  = extno_i + imap - 1 ! 2016-08-16
           obs_npix = getsize_fits(filename, extno = extno_f)       
           ! one partial map (in binary table with explicit pixel indexing)
-          allocate(pixel(0:obs_npix-1), stat = status)
+          allocate(pixel4(0:obs_npix-1), stat = status)
           allocate(signal(0:obs_npix-1), stat = status)
           call read_fits_cut4(filename, int(obs_npix,kind=i4b), &
-               &              pixel, signal, header=hbuffer, units=units1, extno=extno_f)
+               &              pixel4, signal, header=hbuffer, units=units1, extno=extno_f)
           if (present(header) .and. imap == 1) then
              do i=1,nlheader
                 header(i) = hbuffer(i)
              enddo
           endif
           if (present(units)) units(imap) = trim(units1)
-          maxpix = maxval(pixel)
-          minpix = maxval(pixel)
+          maxpix = maxval(pixel4)
+          minpix = maxval(pixel4)
           if (maxpix > (npixtot-1) .or. minpix < 0) then
              print*,'map constructed from file '//trim(filename)//', with pixels in ',minpix,maxpix
              print*,' wont fit in array with ',npixtot,' elements'
              call fatal_error
           endif
           map(:,imap)        = fmiss_effct
-          map(pixel(:),imap) = signal(:)
+          map(pixel4(:),imap) = signal(:)
           imissing = npixtot - obs_npix
           deallocate(signal)
-          deallocate(pixel)
+          deallocate(pixel4)
        enddo
     else 
        print*,'Unable to read the ',nmaps,' required  map(s) from file '//trim(filename)
@@ -2357,3 +2387,578 @@ subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
     return
   end subroutine unfold_weightslist_KLOAD
 
+  !=======================================================================
+  ! read_fits_partial(filename, pixel, cutmap, [header, extno])
+  !   routine to read FITS file with cut sky data : PIXEL, ?, ?, ?, ?, ...
+  !=======================================================================
+#ifndef NO64BITS
+  !=======================================================================
+  subroutine read_fits_partial4_KLOAD(filename, pixel, cutmap, header, extno)
+    !=======================================================================
+    use long_intrinsic, only: long_size
+    character(len=*),                     intent(in)            :: filename
+    integer(I4B),     dimension(0:),      intent(out)           :: pixel
+    real(KMAP),       dimension(0:,1:),   intent(out)           :: cutmap
+    character(len=*), dimension(1:),      intent(out), optional :: header
+    integer(I4B),                         intent(in),  optional :: extno
+    !
+    integer(i8b),     dimension(:), allocatable :: pixel8
+    integer(i8b) :: npix
+    character(len=*), parameter :: code = 'read_fits_partial'
+    !---------------------------------------------------------------------
+    npix = long_size(pixel)
+    if (npix > MAX_I4b) then
+       call fatal_error(code//': PIXEL must be of type I8B to read '//trim(filename))
+    endif
+    allocate(pixel8(0:npix-1))
+    call read_fits_partial8_KLOAD(filename, pixel8, cutmap, header, extno)
+    if (maxval(pixel8) > MAX_I4B) then
+       call fatal_error(code//': PIXEL must be of type I8B to read '//trim(filename))       
+    endif
+    pixel(0:npix-1) = pixel8(0:npix-1)
+    deallocate(pixel8)
+    return
+
+  end subroutine read_fits_partial4_KLOAD
+#endif
+  !=======================================================================
+  subroutine read_fits_partial8_KLOAD(filename, pixel, cutmap, header, extno)
+    !=======================================================================
+    use long_intrinsic, only: long_size
+    character(len=*),                     intent(in)            :: filename
+    integer(I8B),     dimension(0:),      intent(out)           :: pixel
+    real(KMAP),       dimension(0:,1:),   intent(out)           :: cutmap
+    character(len=*), dimension(1:),      intent(out), optional :: header
+    integer(I4B),                         intent(in),  optional :: extno
+
+    integer(I8B) :: obs_npix, npixtot, npix
+
+    integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension
+    integer(I4B) :: blocksize, datacode
+    integer(I4B) :: firstpix, frow, hdutype
+    integer(I4B) :: naxis, nfound, nmove, nrows, nmaps, nmr, i, npix_4
+    integer(I4B) :: readwrite, status, tfields, unit, varidat, width
+    integer(I4B) :: repeat, repeat1, repeat2
+    logical(LGT) :: anynull, extend
+    character(len=20), dimension(MAXDIM) :: ttype, tform, tunit
+    character(len=20) :: extname
+    character(len=80) :: comment
+    real(KMAP) :: nullval
+    !=======================================================================
+
+    npixtot = min(long_size(pixel),long_size(cutmap,1))
+    nmaps   = size(cutmap,2)
+    status=0
+
+    unit = 150
+    nfound = -1
+    anynull = .false.
+
+    readwrite=0
+    call ftopen(unit,filename,readwrite,blocksize,status)
+    if (status > 0) call printerror(status)
+
+    !     determines the presence of image
+    call ftgkyj(unit,'NAXIS', naxis, comment, status)
+    if (status > 0) call printerror(status)
+    if (naxis > 0) then ! there is an image
+       print*,'an image was found in the FITS file '//filename
+       print*,'... it is ignored.'
+    endif
+
+    !     determines the presence of an extension
+    call ftgkyl(unit,'EXTEND', extend, comment, status)
+    if (status > 0) then 
+       print*,'extension expected and not found in FITS file '//filename
+       print*,'abort code'
+       call fatal_error
+    endif
+
+    nmove = +1
+    if (present(extno)) nmove = +1 + extno
+    call ftmrhd(unit, nmove, hdutype, status)
+
+    call assert (hdutype==2, 'this is not a binary table')
+
+    ! reads all the FITS related keywords
+    call ftghbn(unit, MAXDIM, &
+         &        nrows, tfields, ttype, tform, tunit, extname, varidat, &
+         &        status)
+
+    if (.not.  (trim(ttype(1)) == 'PIXEL' &
+         & .or. trim(ttype(1)) == 'PIX'  ) ) call fatal_error('did not find PIXEL column in '//filename)
+    
+
+    if (present(header)) then 
+       header = ""
+       status = 0
+       !call fitstools_mp_get_clean_header(unit, header, filename, status)
+       call get_clean_header(unit, header, filename, status)
+    endif
+
+!     if (present(units)) then 
+!        ! the second column contains the SIGNAL, for which we 
+!        ! already read the units from the header
+!        units = adjustl(tunit(2))
+!     endif
+
+    !        finds the bad data value
+    call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
+    if (status == 202) then ! bad_data not found
+       if (KMAP == SP) nullval = s_bad_value ! default value
+       if (KMAP == DP) nullval = d_bad_value ! default value
+       status = 0
+    endif
+
+    frow = 1
+    firstpix = 1
+    !        parse TFORM keyword to find out the length of the column vector
+    repeat1 = 1
+    repeat2 = 1
+    call ftbnfm(tform(1), datacode, repeat1, width, status)
+    if (tfields > 1) call ftbnfm(tform(2), datacode, repeat2, width, status)
+    repeat = max(repeat1, repeat2)
+
+    !call ftgkyj(unit,'OBS_NPIX',obs_npix,comment,status) ! I4B
+    !call ftgkyk(unit,'OBS_NPIX',obs_npix,comment,status) ! I8B
+    call f90ftgky_(unit,'OBS_NPIX',obs_npix,comment,status) ! generic
+    if (status == 202) then ! obs_npix not found
+       obs_npix = int(nrows, kind=i8b) * repeat
+       status = 0
+    endif
+
+    npix   = min(npixtot, obs_npix)
+    npix_4 = npix
+    nmr    = min(nmaps,   tfields-1)
+    !call ftgcvj   (unit, 1_i4b, frow, firstpix, npix_4, i_bad_value, pixel(0), anynull, status) ! I4B
+    !call ftgcvk   (unit, 1_i4b, frow, firstpix, npix_4, l_bad_value, pixel(0), anynull, status) ! I8B
+    call f90ftgcv_(unit, 1_i4b, frow, firstpix, npix_4, l_bad_value, pixel, anynull, status) ! generic
+    do i=1,nmr
+       call f90ftgcv_(unit, i+1_i4b, frow, firstpix, npix_4, nullval, cutmap(0:npix-1,i), anynull, status)
+    enddo
+
+    !     close the file
+    call ftclos(unit, status)
+
+    !     check for any error, and if so print out error messages
+    if (status > 0) call printerror(status)
+
+    return
+  end subroutine read_fits_partial8_KLOAD
+
+
+  !=======================================================================
+  ! writes a fits file for (polarized) sky data set with columns:
+  !  PIXEL, TEMPERATURE, [Q_POLARISATION, U_POLARISATION]
+  ! 
+  ! write_fits_partial(filename, pixel, cutmap, &
+  !         [, header, coord, nside, order, units, extno])
+  !=======================================================================
+  subroutine write_fits_partial4_KLOAD(filename, pixel, cutmap, &
+       &                     header, coord, nside, order, units, extno)
+    !=======================================================================
+    character(len=*),                   intent(in)           :: filename
+    integer(I4B),     dimension(0:),    intent(in)           :: pixel
+    real(KMAP),       dimension(0:,1:), intent(in)           :: cutmap
+    character(len=*), dimension(1:),    intent(in), optional :: header
+    integer(I4B),                       intent(in), optional :: nside, order
+    character(len=*),                   intent(in), optional :: coord, units
+    integer(I4B),                       intent(in), optional :: extno !, polarisation
+    character(len=*), parameter :: routine = 'write_fits_partial'
+    call sub_write_fits_partial_KLOAD(filename, cutmap, pixel4=pixel, &
+         header=header, coord=coord, nside=nside, order=order, units=units, extno=extno)
+  end subroutine write_fits_partial4_KLOAD
+  !=======================================================================
+#ifndef NO64BITS
+  subroutine write_fits_partial8_KLOAD(filename, pixel, cutmap, &
+       &                     header, coord, nside, order, units, extno)
+    !=======================================================================
+    character(len=*),                   intent(in)           :: filename
+    integer(I8B),     dimension(0:),    intent(in)           :: pixel
+    real(KMAP),       dimension(0:,1:), intent(in)           :: cutmap
+    character(len=*), dimension(1:),    intent(in), optional :: header
+    integer(I4B),                       intent(in), optional :: nside, order
+    character(len=*),                   intent(in), optional :: coord, units
+    integer(I4B),                       intent(in), optional :: extno !, polarisation
+    character(len=*), parameter :: routine = 'write_fits_partial'
+    call sub_write_fits_partial_KLOAD(filename, cutmap, pixel8=pixel, &
+         header=header, coord=coord, nside=nside, order=order, units=units, extno=extno)
+  end subroutine write_fits_partial8_KLOAD
+#endif
+  !=======================================================================
+  subroutine sub_write_fits_partial_KLOAD(filename, cutmap, &
+       &                     pixel4, pixel8, &
+       &                     header, coord, nside, order, units, extno)
+    !=======================================================================
+    use long_intrinsic, only: long_size
+    use pix_tools,      only: nside2npix
+    use misc_utils,     only: assert
+
+    character(len=*),                   intent(in)           :: filename
+    real(KMAP),       dimension(0:,1:), intent(in)           :: cutmap
+    integer(I4B),     dimension(0:),    intent(in), optional :: pixel4
+    integer(I8B),     dimension(0:),    intent(in), optional :: pixel8
+    character(len=*), dimension(1:),    intent(in), optional :: header
+    integer(I4B),                       intent(in), optional :: nside, order
+    character(len=*),                   intent(in), optional :: coord, units
+    integer(I4B),                       intent(in), optional :: extno !, polarisation
+
+    character(len=*), parameter :: routine = 'write_fits_partial'
+    ! --- healpix/fits related variables ----
+    integer(I4B)     :: ncol, grain, nd2
+    integer(I4B)     :: npix_hd, nside_hd, i
+    integer(I4B)     :: maxpix, minpix
+    character(len=1) :: char1, coord_usr
+!    character(len=8) :: char8
+    character(len=20) :: units_usr
+    logical(LGT)     :: done_nside, done_order, done_coord, done_polar, polar_flag
+    integer(I4B)     :: nlheader, extno_i
+
+    ! --- cfitsio related variables ----
+    integer(I4B) ::  status, unit, blocksize, bitpix, naxis, naxes(1)
+    logical(LGT) ::  simple, extend
+    character(LEN=80) :: svalue, comment
+
+    integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension
+    integer(I4B) :: nrows, tfields, varidat
+    integer(I4B) :: frow,  felem, repeat, repeatg, hdutype
+    character(len=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
+    character(len=20), dimension(1:3) :: extnames
+    character(len=80) ::  card
+    character(len=4) :: srepeat, srepeatg
+    character(len=1) :: pform1, pform2
+    character(len=filenamelen) sfilename
+    integer(I8B)     :: obs_npix8, npix_final, obs_npp
+    integer(I4B)     :: obs_npix4, kpix, nside_final
+
+    integer(I4B), save :: nside_old
+    !=======================================================================
+    if (KMAP == SP)  pform2='E'
+    if (KMAP == DP)  pform2='D'
+    if (present(pixel4) .eqv. present(pixel8)) then
+       call fatal_error(routine//': choose either PIXEL4 or PIXEL8')
+    endif
+    if (present(pixel4)) then
+       kpix = I4B
+       pform1='J'
+       obs_npp = long_size(pixel4)
+    elseif (present(pixel8)) then
+       kpix = I8B
+       pform1='K'
+       obs_npp = long_size(pixel8)
+    endif
+    obs_npix8 = long_size(cutmap, 1)
+    nd2       =      size(cutmap, 2)
+    obs_npix4 = int(obs_npix8, kind=I4B)
+    call assert( obs_npix8 == obs_npp, routine//': mismatched size for PIXEL and DATA')
+    ncol = 1 + nd2
+    grain = 1
+    status=0
+    unit = 149
+    blocksize=1
+
+    nlheader = 0
+    if (present(header)) nlheader = size(header)
+    units_usr = ' '
+    if (present(units)) units_usr = units
+    extno_i = 0
+    if (present(extno)) extno_i = extno
+    polar_flag = .false.
+    done_polar = .false.
+!     if (present(polarisation)) then
+!        polar_flag = (polarisation /= 0)
+!        done_polar = .true. ! will ignore the header provided value of POLAR
+!     endif
+
+    if (extno_i == 0) then
+       !*************************************
+       !     create the new empty FITS file
+       !*************************************
+       call ftinit(unit,filename,blocksize,status)
+       if (status > 0) call fatal_error("Error while creating file " &
+            & //trim(filename) &
+            & //". Check path and/or access rights.")
+
+       !     -----------------------------------------------------
+       !     initialize parameters about the FITS image
+       simple=.true.
+       bitpix=32     ! integer*4
+       naxis=0       ! no image
+       naxes(1)=0
+       extend=.true. ! there is an extension
+
+       !     ----------------------
+       !     primary header
+       !     ----------------------
+       !     write the required header keywords
+       call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
+
+       !     write the current date
+       call ftpdat(unit,status) ! format ccyy-mm-dd
+
+    else
+       !*********************************************
+       !     reopen an existing file and go to the end
+       !*********************************************
+       ! remove the leading '!' (if any) when reopening the same file
+       sfilename = adjustl(filename)
+       if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
+       call ftopen(unit,sfilename,1_i4b,blocksize,status)
+       ! move to last extension written
+       call ftmahd(unit, 1_i4b+ extno_i, hdutype, status)
+
+    endif
+
+       
+    !     ----------------------
+    !     new extension
+    !     ----------------------
+    call ftcrhd(unit, status)
+
+       !     writes required keywords
+       !     repeat = 1024
+    repeat = 1
+    if (obs_npix8 < repeat) repeat = 1
+    nrows    = (obs_npix8 + repeat - 1)/ repeat ! naxis1
+    write(srepeat,'(i4)') repeat
+    srepeat = adjustl(srepeat)
+
+    repeatg = repeat * grain
+    write(srepeatg,'(i4)') repeatg
+    srepeatg = adjustl(srepeatg)
+
+    tfields  = ncol
+    ttype(1) = 'PIXEL'
+    if (nd2  == 1) then
+       done_polar = .true.
+       polar_flag = .false.
+       ttype(2)   = 'TEMPERATURE  '
+    elseif (nd2  == 3) then
+       done_polar = .true.
+       polar_flag = .true.
+       ttype(2:4) = (/'TEMPERATURE   ' , & 
+            &         'Q_POLARISATION' , &
+            &         'U_POLARISATION' /)
+    else
+       do i=2,ncol
+          ttype(i) = 'C'//string(i-1, format='(i2.2)') ! C01, C02, ...
+       enddo
+    endif
+
+    tform(1)      = trim(srepeat) //pform1
+    tform(2:ncol) = trim(srepeatg)//pform2
+
+    tunit =  ' '      ! optional, will not appear
+    tunit(2:ncol) = units_usr
+
+    extname  = 'SKY_OBSERVATION'      ! default, will be overridden by user provided one if any
+    ! if (polar_flag) extname = extnames(1+extno_i)
+    varidat  = 0
+    call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
+         &     extname, varidat, status)
+
+    call ftpcom(unit,'------------------------------------------',status)
+    call ftpcom(unit,'          Pixelisation Specific Keywords    ',status)
+    call ftpcom(unit,'------------------------------------------',status)
+    call ftpkys(unit,'PIXTYPE','HEALPIX ',' HEALPIX Pixelisation',status)      
+    call ftpkyu(unit,'NSIDE',   ' ',status) ! place holder, will be updated later on
+    call ftpkyu(unit,'ORDERING',' ',status)
+    call ftpkys(unit,'COORDSYS',' ',' ',status)
+    call ftpcom(unit,'  G = Galactic, E = ecliptic, C = celestial = equatorial', status)   
+    call ftpkyl(unit,'POLAR',polar_flag,'Polarisation included in file (T/F)',status)
+    call ftpcom(unit,'------------------------------------------',status)
+    call ftpcom(unit,'          Data Specific Keywords    ',status)
+    call ftpcom(unit,'------------------------------------------',status)
+    call ftpkys(unit,'INDXSCHM','EXPLICIT',' Indexing : IMPLICIT or EXPLICIT', status)
+!     call ftpkyj(unit,'GRAIN',  grain,     ' Grain of pixel indexing',status)
+!     call ftpcom(unit,'GRAIN=0 : no indexing of pixel data                         (IMPLICIT)',status)
+!     call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data                     (EXPLICIT)',status)
+!     call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
+    call ftpkys(unit,'OBJECT','PARTIAL ',' Sky coverage represented by data',status)   
+    if (KPIX == I4B) then
+       call ftpkyj(unit,'OBS_NPIX',obs_npix4, ' Number of pixels observed and recorded',status)
+    else
+       call ftpkyk(unit,'OBS_NPIX',obs_npix8, ' Number of pixels observed and recorded',status)
+    endif
+
+    ! add required Healpix keywords (NSIDE, ORDER) if provided by user
+    done_order = .false.
+    if (present(order)) then
+       !        char8 = order
+       !        call ftupch(char8)
+       !        if (char8(1:4) == 'RING') then
+       if (order == 1) then
+          call ftukys(unit, 'ORDERING','RING',' Pixel ordering scheme, either RING or NESTED',status)
+          done_order = .true.
+          !        elseif (char8(1:4) == 'NEST') then
+       elseif (order == 2) then
+          call ftukys(unit, 'ORDERING','NESTED',' Pixel ordering scheme, either RING or NESTED ',status)
+          done_order = .true.
+       else
+          print*,'Invalid ORDER given : ',order, ' instead of 1 (RING) or 2 (NESTED)'
+       endif
+    endif
+
+    done_nside = .false.
+    if (present(nside)) then
+       if (nside2npix(nside) > 0) then ! valid nside
+          call ftukyj(unit,'NSIDE',nside,' Resolution parameter for HEALPIX', status)
+          done_nside = .true.
+          nside_final = nside
+       else
+          print*,'Invalid NSIDE given : ',nside
+       endif
+    endif
+
+    ! add non required Healpix keyword (COORD)
+    done_coord = .false.
+    if (present(coord)) then
+       coord_usr = adjustl(coord)
+       char1 = coord_usr(1:1)
+       call ftupch(char1) ! uppercase
+       if (char1 == 'C' .or. char1 == 'Q') then
+          coord_usr = 'C'
+          done_coord = .true.
+       elseif (char1 == 'G') then
+          coord_usr = 'G'
+          done_coord = .true.
+       elseif (char1 == 'E' ) then
+          coord_usr='E'
+          done_coord = .true.
+       else
+          print*,'Unrecognised COORD given : ',coord,' instead of C, G, or E'
+          print*,'Proceed at you own risks '
+          coord_usr = char1
+          done_coord = .true.
+       endif
+       if (done_coord) then
+          call ftukys(unit, 'COORDSYS',coord_usr,' Pixelisation coordinate system',status)
+       endif
+    endif
+
+
+    !    write the user provided header literally, except for  PIXTYPE, TFORM*, TTYPE*, TUNIT*, INDXSCHM and GRAIN
+    !    copy NSIDE, ORDERING and COORDSYS and POLAR if they are valid and not already given
+    do i=1,nlheader
+       card = header(i)
+       if (card(1:5) == 'TTYPE' .or. card(1:5) == 'TFORM' .or. card(1:7) == 'PIXTYPE') then
+          continue ! don't keep them
+       else if (card(1:8) == 'INDXSCHM') then
+          continue
+       else if (card(1:5) == 'GRAIN') then ! already written above
+          continue
+       else if (card(1:13) == 'COMMENT GRAIN' .or. card(1:14) == 'COMMENT  GRAIN') then ! already written above
+          continue
+       else if (card(1:5) == 'TUNIT') then 
+          if (trim(units_usr) == '') then
+             call ftucrd(unit,'TUNIT'//card(6:7),card, status) !update TUNIT2 and above
+          endif
+       else if (card(1:5) == 'NSIDE') then
+          call ftpsvc(card, svalue, comment, status)
+          read(svalue,*) nside_hd
+          npix_hd = nside2npix(nside_hd)
+          if (.not. done_nside .and. npix_hd > 0) then
+             call ftucrd(unit,'NSIDE',card, status) !update NSIDE
+             done_nside = .true.
+             nside_final = nside_hd
+          endif
+       else if (card(1:8) == 'ORDERING') then
+          call ftpsvc(card, svalue, comment, status)
+          svalue = adjustl(svalue)
+          svalue = svalue(2:8) ! remove leading quote
+          call ftupch(svalue)
+          if (.not. done_order .and. (svalue(1:4)=='RING' .or. svalue(1:6) == 'NESTED')) then
+             call ftucrd(unit,'ORDERING',card, status) !update ORDERING
+             done_order = .true.
+          endif
+       else if (card(1:8) == 'COORDSYS') then
+          if (.not. done_coord) call putrec(unit,card, status)
+          done_coord = .true.
+       else if (card(1:5) == 'POLAR') then
+          if (.not. done_polar) then
+             call ftucrd(unit,'POLAR', card, status) ! update POLAR
+             done_polar = .true.
+          endif
+       else if (card(1:7) == 'EXTNAME') then
+          call ftucrd(unit, 'EXTNAME', card, status)
+       else if (card(1:1) /= ' ') then ! if none of the above, copy to FITS file
+          call putrec(unit, card, status)
+       endif
+10     continue
+    enddo
+
+
+    ! test that required keywords have been provided in some way
+    if (.not. done_nside) then
+       print*,routine//'> NSIDE is a Healpix required keyword, '
+       print*,routine//'>  it was NOT provided either as routine argument or in the input header'
+       print*,routine//'>  abort execution, file not written'
+       call fatal_error
+    endif
+    if (.not. done_order) then
+       print*,routine//'> ORDER is a Healpix required keyword, '
+       print*,routine//'>  it was NOT provided either as routine argument or in the input header'
+       print*,routine//'>  abort execution, file not written'
+       call fatal_error
+    endif
+!     if ((.not. done_polar) .and. extno_i >=2) then
+!        print*,routine//'> Warning: POLAR keyword not set while 3 extensions have been written'
+!     endif
+
+    ! check that NSIDE is the same for all extensions
+    if (extno_i == 0) then
+       nside_old = nside_final
+    else
+       if (nside_final /= nside_old) then
+          print*,routine//'> Inconsistent NSIDE: ',nside_final, nside_old
+          print*,routine//'> Should use same NSIDE for all extensions'
+       endif
+    endif
+
+    ! check validity of PIXEL
+    npix_final = nside2npix(nside_final)
+    if (kpix == I4B) then
+       minpix = minval(pixel4(0:obs_npix8-1))
+       maxpix = maxval(pixel4(0:obs_npix8-1))
+    else
+       minpix = minval(pixel8(0:obs_npix8-1))
+       maxpix = maxval(pixel8(0:obs_npix8-1))
+    endif
+    if (minpix < 0 .or. maxpix > npix_final-1) then
+       print*,routine//'> Actual pixel range = ',minpix,maxpix
+       print*,routine//'> expected range (for Nside =',nside_final,') : 0, ',npix_final-1
+       print*,routine//'> ABORT execution, file not written.'
+       call fatal_error
+    endif
+    if (obs_npix8 > npix_final) then 
+       print*,routine//'> The actual number of pixels ',obs_npix8
+       print*,routine//'> is larger than the number of pixels over the whole sky : ',npix_final
+       print*,routine//'> for Nside = ',nside_final
+       print*,routine//'> ABORT execution, file not written.'
+       call fatal_error
+    endif
+
+    !     write the extension one column by one column
+    frow   = 1  ! starting position (row)
+    felem  = 1  ! starting position (element)
+    if (kpix == I4B) then
+       call f90ftpcl_(unit, 1_i4b, frow, felem, obs_npix4, pixel4, status)
+    else
+       call f90ftpcl_(unit, 1_i4b, frow, felem, obs_npix4, pixel8, status)
+    endif
+    do i=1, nd2
+       call f90ftpcl_(unit, i+1_i4b, frow, felem, obs_npix4, cutmap(0:,i), status)
+    enddo
+
+    !     ----------------------
+    !     close and exit
+    !     ----------------------
+
+    !     close the file and free the unit number
+    call ftclos(unit, status)
+
+    !     check for any error, and if so print out error messages
+    if (status > 0) call printerror(status)
+
+    return
+  end subroutine sub_write_fits_partial_KLOAD
diff --git a/src/f90/mod/fitstools.F90 b/src/f90/mod/fitstools.F90
index 2f3d3ce..6705185 100644
--- a/src/f90/mod/fitstools.F90
+++ b/src/f90/mod/fitstools.F90
@@ -84,6 +84,8 @@ module fitstools
   ! subroutine read_bintod
   ! subroutine write_bintabh
   ! subroutine unfold_weights
+  ! subroutine read_fits_partial
+  ! subroutine write_fits_partial
   ! ----------------------------------
   !
   ! subroutine read_fits_cut4               ?
@@ -119,6 +121,7 @@ module fitstools
   real(kind=SP),     private, parameter :: s_bad_value = HPX_SBADVAL
   real(kind=DP),     private, parameter :: d_bad_value = HPX_DBADVAL
   integer(kind=I4B), private, parameter :: i_bad_value = -1637500000
+  integer(kind=I8B), private, parameter :: l_bad_value = -1637500000_I8B
   integer(I4B) ,     private, parameter :: nchunk_max  = 12000
   integer(I4B),      private, parameter :: MAXDIM_TOP  = 199 ! < 999
 
@@ -167,6 +170,24 @@ module fitstools
 #endif
   end interface
 
+  interface read_fits_partial
+#ifdef NO64BITS
+     module procedure read_fits_partial8_s, read_fits_partial8_d
+#else
+     module procedure read_fits_partial8_s, read_fits_partial8_d, &
+          & read_fits_partial4_s, read_fits_partial4_d
+#endif
+  end interface
+
+  interface write_fits_partial
+#ifdef NO64BITS
+     module procedure write_fits_partial4_s, write_fits_partial4_d
+#else
+     module procedure write_fits_partial8_s, write_fits_partial8_d, &
+          & write_fits_partial4_s, write_fits_partial4_d
+#endif
+  end interface
+
   !
   interface input_tod
      module procedure input_tod_s,input_tod_d
@@ -217,11 +238,19 @@ module fitstools
   end interface
 
   interface f90ftpcl_
-     module procedure f90ftpcle, f90ftpcld
+#ifdef NO64BITS
+     module procedure f90ftpcle, f90ftpcld, f90ftpclj
+#else
+     module procedure f90ftpcle, f90ftpcld, f90ftpclj, f90ftpclk
+#endif
   end interface
 
   interface f90ftgcv_
-     module procedure f90ftgcve, f90ftgcvd
+#ifdef NO64BITS
+     module procedure f90ftgcve, f90ftgcvd, f90ftgcvj
+#else
+     module procedure f90ftgcve, f90ftgcvd, f90ftgcvj, f90ftgcvk
+#endif
   end interface
 
   interface f90ftgpv_
@@ -229,7 +258,11 @@ module fitstools
   end interface
 
   interface f90ftgky_
-     module procedure f90ftgkye, f90ftgkyd
+#ifdef NO64BITS
+     module procedure f90ftgkye, f90ftgkyd, f90ftgkyj
+#else
+     module procedure f90ftgkye, f90ftgkyd, f90ftgkyj, f90ftgkyk
+#endif
   end interface
 
   interface map_bad_pixels
@@ -249,6 +282,7 @@ module fitstools
   public :: read_fits_cut4, write_fits_cut4, & 
        & input_map, read_bintab,  &
        & output_map, write_bintab
+  public :: read_fits_partial, write_fits_partial
   public :: write_plm
   public :: fits2cl, read_asctab, write_asctab
   public :: read_dbintab, write_dbintab
@@ -262,7 +296,7 @@ module fitstools
 contains
 
   !-------------------------------------------------------------------------------
-  ! generic interface F90FTPCL_ for FITSIO's FTPCLE and FTPCLD
+  ! generic interface F90FTPCL_ for FITSIO's FTPCL[E,D,J,K]
   !           writes data in ASCTAB or BINTAB
   subroutine f90ftpcle(unit, colnum, frow, felem, np, data, status)
     integer(I4B), intent(in)  :: unit, colnum, frow, felem, np
@@ -278,8 +312,24 @@ contains
     call ftpcld(unit, colnum, frow, felem, np, data, status)
     return
   end subroutine f90ftpcld
+  subroutine f90ftpclj(unit, colnum, frow, felem, np, data, status)
+    integer(I4B), intent(in)  :: unit, colnum, frow, felem, np
+    integer(I4B), intent(out) :: status
+    integer(I4B), intent(in), dimension(0:)  :: data
+    call ftpclj(unit, colnum, frow, felem, np, data, status)
+    return
+  end subroutine f90ftpclj
+#ifndef NO64BITS
+  subroutine f90ftpclk(unit, colnum, frow, felem, np, data, status)
+    integer(I4B), intent(in)  :: unit, colnum, frow, felem, np
+    integer(I4B), intent(out) :: status
+    integer(I8B), intent(in), dimension(0:)  :: data
+    call ftpclk(unit, colnum, frow, felem, np, data, status)
+    return
+  end subroutine f90ftpclk
+#endif
   !-------------------------------------------------------------------------------
-  ! generic interface F90FTGCV_ for FITSIO's FTGCVE and FTGCVD
+  ! generic interface F90FTGCV_ for FITSIO's FTGCV[E,D,J,K]
   !           reads data from BINTAB
   subroutine f90ftgcve(unit, colnum, frow, felem, np, nullval, data, anynull, status)
     integer(I4B), intent(in)  :: unit, colnum, frow, felem, np
@@ -299,6 +349,26 @@ contains
     call ftgcvd(unit, colnum, frow, felem, np, nullval, data, anynull, status)
     return
   end subroutine f90ftgcvd
+  subroutine f90ftgcvj(unit, colnum, frow, felem, np, nullval, data, anynull, status)
+    integer(I4B), intent(in)  :: unit, colnum, frow, felem, np
+    integer(I4B), intent(out) :: status
+    logical(LGT), intent(out) :: anynull
+    integer(I4B), intent(out), dimension(0:) :: data
+    integer(I4B), intent(in)                 :: nullval
+    call ftgcvj(unit, colnum, frow, felem, np, nullval, data, anynull, status)
+    return
+  end subroutine f90ftgcvj
+#ifndef NO64BITS
+  subroutine f90ftgcvk(unit, colnum, frow, felem, np, nullval, data, anynull, status)
+    integer(I4B), intent(in)  :: unit, colnum, frow, felem, np
+    integer(I4B), intent(out) :: status
+    logical(LGT), intent(out) :: anynull
+    integer(I8B), intent(out), dimension(0:) :: data
+    integer(I8B), intent(in)                 :: nullval
+    call ftgcvk(unit, colnum, frow, felem, np, nullval, data, anynull, status)
+    return
+  end subroutine f90ftgcvk
+#endif
   !-------------------------------------------------------------------------------
   ! generic interface F90FTGPV_ for FITSIO's FTGPVE and FTGPVD
   !           reads data from IMAGE
@@ -321,7 +391,7 @@ contains
     return
   end subroutine f90ftgpvd
   !-------------------------------------------------------------------------------
-  ! generic interface F90FTGKY_ for FITSIO's FTGKYE and FTGKYD
+  ! generic interface F90FTGKY_ for FITSIO's FTGKY[E,D,J,K]
   !           reads a keyword
   subroutine f90ftgkye(unit, keyword, value, comment, status)
     integer(I4B),     intent(in)  :: unit
@@ -341,14 +411,34 @@ contains
     call ftgkyd(unit, keyword, value, comment, status)
     return
   end subroutine f90ftgkyd
+  subroutine f90ftgkyj(unit, keyword, value, comment, status)
+    integer(I4B),     intent(in)  :: unit
+    character(len=*), intent(in)  :: keyword
+    integer(I4B),     intent(out) :: status
+    character(len=*), intent(out) :: comment
+    integer(I4B),     intent(out) :: value
+    call ftgkyj(unit, keyword, value, comment, status)
+    return
+  end subroutine f90ftgkyj
+#ifndef NO64BITS
+  subroutine f90ftgkyk(unit, keyword, value, comment, status)
+    integer(I4B),     intent(in)  :: unit
+    character(len=*), intent(in)  :: keyword
+    integer(I4B),     intent(out) :: status
+    character(len=*), intent(out) :: comment
+    integer(I8B),     intent(out) :: value
+    call ftgkyk(unit, keyword, value, comment, status)
+    return
+  end subroutine f90ftgkyk
+#endif
   !-------------------------------------------------------------------------------
 
 
   ! define routine with SP I/O
-#include "fits_s_inc.f90"
+#include "fits_s_inc.F90"
 
   ! define routine with DP I/O
-#include "fits_d_inc.f90"
+#include "fits_d_inc.F90"
 
 
 
@@ -402,14 +492,14 @@ contains
     call ftgkyj(unit,'NAXIS', naxis, comment, status)
     if (status > 0) call printerror(status)
     if (naxis > 0) then ! there is an image
-       print*,'an image was found in the FITS file '//filename
+       print*,'an image was found in the FITS file '//trim(filename)
        print*,'... it is ignored.'
     endif
 
     !     determines the presence of an extension
     call ftgkyl(unit,'EXTEND', extend, comment, status)
     if (status > 0) then 
-       print*,'extension expected and not found in FITS file '//filename
+       print*,'extension expected and not found in FITS file '//trim(filename)
        print*,'abort code'
        call fatal_error
     endif
@@ -425,7 +515,7 @@ contains
          &        nrows, tfields, ttype, tform, tunit, extname, varidat, &
          &        status)
     if (tfields < 4) then
-       print*,'Expected 4 columns in FITS file '//filename
+       print*,'Expected 4 columns in FITS file '//trim(filename)
        print*,'found ',tfields
        if (tfields < 2) call fatal_error
        if (.not.  (trim(ttype(1)) == 'PIXEL' &
@@ -614,8 +704,8 @@ contains
        !     writes required keywords
        !     repeat = 1024
     repeat = 1
-    nrows    = (obs_npix + repeat - 1)/ repeat ! naxis1
     if (obs_npix < repeat) repeat = 1
+    nrows    = (obs_npix + repeat - 1)/ repeat ! naxis1
     write(srepeat,'(i4)') repeat
     srepeat = adjustl(srepeat)
 
@@ -633,7 +723,7 @@ contains
     tunit =  ' '      ! optional, will not appear
     tunit(2) = units_usr
     tunit(4) = units_usr
-    extname  = 'SKY_OBSERVATION'      ! default, will be overide by user provided one if any
+    extname  = 'SKY_OBSERVATION'      ! default, will be overridden by user provided one if any
     if (polar_flag) extname = extnames(1+extno_i)
     varidat  = 0
     call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
@@ -652,10 +742,10 @@ contains
     call ftpcom(unit,'          Data Specific Keywords    ',status)
     call ftpcom(unit,'------------------------------------------',status)
     call ftpkys(unit,'INDXSCHM','EXPLICIT',' Indexing : IMPLICIT or EXPLICIT', status)
-    call ftpkyj(unit,'GRAIN',  grain,     ' Grain of pixel indexing',status)
-    call ftpcom(unit,'GRAIN=0 : no indexing of pixel data (IMPLICIT) ',status)
-    call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)',status)
-    call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
+!     call ftpkyj(unit,'GRAIN',  grain,     ' Grain of pixel indexing',status)
+!     call ftpcom(unit,'GRAIN=0 : no indexing of pixel data                         (IMPLICIT)',status)
+!     call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data                     (EXPLICIT)',status)
+!     call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
     call ftpkys(unit,'OBJECT','PARTIAL ',' Sky coverage represented by data',status)     
     call ftpkyj(unit,'OBS_NPIX',obs_npix, ' Number of pixels observed and recorded',status)
 
@@ -716,7 +806,7 @@ contains
     endif
 
 
-    !    write the user provided header literally, except for  PIXTYPE, TFORM*, TTYPE*, TUNIT* and INDXSCHM
+    !    write the user provided header literally, except for  PIXTYPE, TFORM*, TTYPE*, TUNIT*, INDXSCHM and GRAIN
     !    copy NSIDE, ORDERING and COORDSYS and POLAR if they are valid and not already given
     do i=1,nlheader
        card = header(i)
@@ -724,6 +814,10 @@ contains
           continue ! don't keep them
        else if (card(1:8) == 'INDXSCHM') then
           continue
+       else if (card(1:5) == 'GRAIN') then ! already written above
+          continue
+       else if (card(1:13) == 'COMMENT GRAIN' .or. card(1:14) == 'COMMENT  GRAIN') then ! already written above
+          continue
        else if (card(1:5) == 'TUNIT') then 
           if ((card(6:6) == '2' .or. card(6:6) == '4') .and. trim(units_usr) == '') then
              call ftucrd(unit,'TUNIT'//card(6:6),card, status) !update TUNIT2 and TUNIT4
@@ -1968,7 +2062,8 @@ contains
     INTEGER(I4B) :: status,unit,readwrite,blocksize,naxis
     CHARACTER(LEN=80) :: comment, ttype1
     LOGICAL(LGT) ::  extend, anyf
-    INTEGER(I4B)::  nmove, hdutype, idmax, nrows
+    INTEGER(I4B)::  nmove, hdutype, nrows
+    integer(I4B), dimension(1) :: idmax ! because ftgcvj expects an array (2020-08-24)
 
     !-----------------------------------------------------------------------
     status=0
@@ -2015,9 +2110,9 @@ contains
           ttype1 = trim(strupcase(adjustl(ttype1)))
           if (trim(ttype1(1:5)) == 'INDEX') then
              call ftgkyj(unit, 'NAXIS2', nrows, comment, status) ! find number of rows
-             call ftgcvj(unit, 1_i4b, nrows, 1_i4b, 1_i4b, 0_i4b, idmax, anyf, status) ! read element on last row of first column
+             call ftgcvj(unit, 1_i4b, nrows, 1_i4b, 1_i4b, 0_i4b, idmax(1), anyf, status) ! read element on last row of first column
              if (status == 0) then
-                lmax = int(sqrt(   real(idmax-1, kind = DP)  ) )
+                lmax = int(sqrt(   real(idmax(1)-1, kind = DP)  ) )
                 if (lmax > 0) goto 1000
              endif
           endif
@@ -2759,28 +2854,29 @@ contains
     character(len=*), parameter :: primer_url = 'http://healpix.sf.net/pdf/intro.pdf'
     !====================================================================
     
+    !n_ext = getnumext_fits(mapfile)
     npixtot = getsize_fits(mapfile, nmaps = nmaps, ordering=order_map, nside=nsmax,&    
          &              mlpol=mlpol, type = type, polarisation = polar_fits, &
          &             coordsys=coordsys, polcconv=polcconv)
     
     if (nsmax<=0) then
-       print*,"Keyword NSIDE not found in FITS header!"
+       print*,"Keyword NSIDE not found in FITS header of "//trim(mapfile)//" !"
        call fatal_error(code)
     endif
     if (type == 3) npixtot = nside2npix(nsmax) ! cut sky input data set
     if (nsmax/=npix2nside(npixtot)) then
        print 9000,"FITS header keyword NSIDE does not correspond"
-       print 9000,"to the size of the map!"
+       print 9000,"to the size of the map in "//trim(mapfile)//" !"
        call fatal_error(code)
     endif
 
     if (polarisation .and. (nmaps >=3) .and. polar_fits == -1) then
-       print 9000,"The input fits file MAY NOT contain polarisation data."
+       print 9000,"The input FITS file "//trim(mapfile)//" MAY NOT contain polarisation data."
        print 9000,"Proceed at your own risk"
     endif
     
     if (polarisation .and. (nmaps<3 .or. polar_fits ==0)) then
-       print 9000,"The file does NOT contain polarisation maps"
+       print 9000,"The FITS file "//trim(mapfile)//" does NOT contain polarisation maps"
        print 9000,"only the temperature field will be analyzed"
        polarisation = .false.
     endif
@@ -2809,7 +2905,7 @@ contains
     !     --- check ordering scheme ---
     if ((order_map/=1).and.(order_map/=2)) then
        print 9000,"The ordering scheme of the map must be RING or NESTED."
-       print 9000,"No ordering specification is given in the FITS-header!"
+       print 9000,"No ordering specification is given in the FITS-header of "//trim(mapfile)//" !"
        call fatal_error(code)
     endif
     
diff --git a/src/f90/mod/gen_fits_code b/src/f90/mod/gen_fits_code
index a17ab38..bc5fed4 100755
--- a/src/f90/mod/gen_fits_code
+++ b/src/f90/mod/gen_fits_code
@@ -3,9 +3,9 @@
 # EH, Dec 2004 -- Jan 2005
 
 SetDef(){
-template_file=fits_template.f90
-outfile_s=fits_s_inc.f90
-outfile_d=fits_d_inc.f90
+template_file=fits_template.F90
+outfile_s=fits_s_inc.F90
+outfile_d=fits_d_inc.F90
 }
 
 #--------
diff --git a/src/f90/mod/head_fits.F90 b/src/f90/mod/head_fits.F90
index 477ff6f..b804850 100644
--- a/src/f90/mod/head_fits.F90
+++ b/src/f90/mod/head_fits.F90
@@ -524,7 +524,7 @@ contains
           endif
        endif
     enddo
-    if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+    if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
     if (present(comment)) comment = ' '
     if (present(count))   count = count_in
     return
@@ -563,7 +563,7 @@ contains
           endif
        endif
     enddo
-    if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+    if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
     if (present(comment)) comment = ' '
     if (present(count))   count = count_in
     return
@@ -673,7 +673,7 @@ contains
           return
        endif
     enddo
-    if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+    if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
     if (present(comment)) comment = ' '
     if (present(count))   count = count_in
     return
@@ -712,7 +712,7 @@ contains
           endif
        endif
     enddo
-    if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+    if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
     if (present(comment)) comment = ' '
     if (present(count))   count = count_in
     return
@@ -753,7 +753,7 @@ contains
           endif
        endif
     enddo
-    if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+    if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
     if (present(comment)) comment = ' '
     if (present(count))   count = count_in
     return
@@ -792,7 +792,7 @@ contains
           endif
        endif
     enddo
-    if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+    if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
     if (present(comment)) comment = ' '
     if (present(count))   count = count_in
     return
@@ -878,14 +878,15 @@ contains
     character(len=80), dimension(1:20) :: my_units
     character(len=8) :: my_ordering
     character(len=FILENAMELEN) :: udtype, my_beam_leg
-    logical(LGT) :: do_full, do_alm, do_cl, do_cut, do_polar, &
-         & do_bcross, do_reset, do_asym
+    logical(LGT) :: do_full, do_alm, do_cl, do_cut, do_partial, &
+         & do_polar, do_bcross, do_reset, do_asym, do_map
 
     character(len=*), parameter :: code = 'write_minimal_header'
     !--------------------------------------------------------------------------------
     ! find out type of data to be contained in FITS file
     udtype = strupcase(dtype)
-    do_alm = .false. ; do_cl  = .false. ; do_full = .false. ; do_cut = .false.
+    do_alm = .false. ; do_cl  = .false. ; do_full = .false. 
+    do_cut = .false. ; do_partial = .false.
     select case (trim(udtype))
     case ('ALM')
        do_alm = .true.
@@ -895,11 +896,14 @@ contains
        do_full = .true.
     case ('CUTMAP')
        do_cut = .true.
+    ! case ('PARTIAL')
+       do_partial = .true.
     case default
        print*,'Invalid choice of datatype:'//trim(dtype)
        print*,'Should be one of: ALM, alm, CL, cl, CUTMAP, cutmap, MAP, map'
        call fatal_error(code)
     end select
+    do_map = (do_cut .or. do_full .or. do_partial)
 
     !
     if (present(order) .and. present(ordering)) then
@@ -908,7 +912,7 @@ contains
     endif
 
     ! check for required keywords for map
-    if (do_cut .or. do_full) then
+    if (do_map) then
        if (.not. present(nside)) then
           print*,'NSIDE required for maps in '//code
           call fatal_error(code)
@@ -1108,7 +1112,7 @@ contains
 
     endif ! cl
 
-    if (do_full .or. do_cut) then
+    if (do_map) then
        !--------------------------------------------------------------------------------------------
        !                            minimal header for FULL SKY MAP and CUT SKY MAP
        !--------------------------------------------------------------------------------------------
@@ -1156,27 +1160,54 @@ contains
           call add_card(header,"COMMENT","Full sky data")
           call add_card(header,"OBJECT","FULLSKY")
           call add_card(header,"INDXSCHM","IMPLICIT"," Indexing : IMPLICIT or EXPLICIT")
-          call add_card(header,"GRAIN", 0, " Grain of pixel indexing")
+!          call add_card(header,"GRAIN", 0, " Grain of pixel indexing")
        endif
-       if (do_cut)  then
+       if (do_cut .or. do_partial)  then
           call add_card(header,"EXTNAME","'CUT SKY MAP'",update=.true.)
           call add_card(header,"COMMENT","Cut sky data")
           call add_card(header,"OBJECT","PARTIAL")
           call add_card(header,"INDXSCHM","EXPLICIT"," Indexing : IMPLICIT or EXPLICIT")
-          call add_card(header,"GRAIN", 1, " Grain of pixel indexing")
+!          call add_card(header,"GRAIN", 1, " Grain of pixel indexing")
        endif
-       call add_card(header,"COMMENT","GRAIN=0 : no indexing of pixel data                         (IMPLICIT)")
-       call add_card(header,"COMMENT","GRAIN=1 : 1 pixel index -> 1 pixel data                     (EXPLICIT)")
-       call add_card(header,"COMMENT","GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)")
+!        call add_card(header,"COMMENT","GRAIN=0 : no indexing of pixel data                         (IMPLICIT)")
+!        call add_card(header,"COMMENT","GRAIN=1 : 1 pixel index -> 1 pixel data                     (EXPLICIT)")
+!        call add_card(header,"COMMENT","GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)")
        call add_card(header) ! blank line
        call add_card(header,"POLAR",do_polar," Polarisation included (True/False)")
        if (do_cut) then
+          call add_card(header) ! blank line
+          call add_card(header,"TTYPE1", "PIXEL","pixel index")
+          call add_card(header,"TUNIT1", "",     "")
+          call add_card(header)
           if (present(units)) then
              my_units = units
-             call add_card(header,"TUNIT2", my_units(1),"physical unit of signal map")
-             call add_card(header,"TUNIT4", my_units(1),"physical unit of error map")
+             call add_card(header,"TUNIT2", my_units(1),"physical units of map")
+             call add_card(header,"TUNIT4", my_units(1),"physical units of map")
           endif
        endif
+!        if (do_partial) then
+!           if (present(units)) my_units = units
+!           it = 1 ; iq = 1 ; iu = 1
+!           call add_card(header) ! blank line
+!           call add_card(header,"TTYPE1", "PIXEL","pixel index")
+!           call add_card(header,"TUNIT1", "",     "")
+!           call add_card(header)
+
+!           call add_card(header) ! blank line
+!           call add_card(header,"TTYPE2", "TEMPERATURE","Temperature map")
+!           call add_card(header,"TUNIT2", my_units(it), "physical units of map")
+!           call add_card(header)
+
+!           if (do_polar) then
+!              call add_card(header,"TTYPE3", "Q_POLARISATION","Q Polarisation map")
+!              call add_card(header,"TUNIT3", my_units(iq),    "physical units of map")
+!              call add_card(header)
+
+!              call add_card(header,"TTYPE4", "U_POLARISATION","U Polarisation map")
+!              call add_card(header,"TUNIT4", my_units(iu),    "physical units unit")
+!              call add_card(header)
+!           endif
+!        endif
        if (do_full) then
           call add_card(header,"DERIV",my_deriv," Derivative included (0, 1 or 2)")
 
@@ -1199,12 +1230,10 @@ contains
           call add_card(header)
 
           if (do_polar) then
-             !call add_card(header,"TTYPE2", "Q-POLARISATION","Q Polarisation map")
              call add_card(header,"TTYPE2", "Q_POLARISATION","Q Polarisation map")
              call add_card(header,"TUNIT2", my_units(iq),"map unit")
              call add_card(header)
 
-             !call add_card(header,"TTYPE3", "U-POLARISATION","U Polarisation map")
              call add_card(header,"TTYPE3", "U_POLARISATION","U Polarisation map")
              call add_card(header,"TUNIT3", my_units(iu),"map unit")
              call add_card(header)
diff --git a/src/f90/mod/healpix_types.F90 b/src/f90/mod/healpix_types.F90
index fc4dff8..a4fe0e6 100644
--- a/src/f90/mod/healpix_types.F90
+++ b/src/f90/mod/healpix_types.F90
@@ -39,7 +39,7 @@ MODULE healpix_types
   !            Mars 2008: i8b same as i4b on machines not supporting 64 bits (NO64BITS flag set)
   !            Feb  2009: introduce healpix_version
   !
-  character(len=*), PARAMETER, public :: healpix_version = '3.60'
+  character(len=*), PARAMETER, public :: healpix_version = '3.81'
   INTEGER, PARAMETER, public :: i4b = SELECTED_INT_KIND(9)
 #ifdef NO64BITS
   INTEGER, PARAMETER, public :: i8b = i4b
diff --git a/src/f90/mod/pix_tools.F90 b/src/f90/mod/pix_tools.F90
index f9ec2d0..0df526a 100644
--- a/src/f90/mod/pix_tools.F90
+++ b/src/f90/mod/pix_tools.F90
@@ -871,9 +871,11 @@ contains
     integer(i4b), dimension(1:,1:), intent(inout) :: ringphi
     integer(i4b),                   intent(inout) :: ngr
 
-    integer(i4b) :: i, j, k, kk, ngr_out, diff, iphi, i0, nrh
+    !integer(i4b) :: i, j, k, kk, ngr_out, diff, iphi, i0, nrh
+    integer(i4b) :: i, j, k, kk, ngr_out, iphi, i0, nrh, nshrink, np
     real(dp), dimension(1:2*nsboost+1) :: phiw, phie
     real(dp) :: dd, dph, phic
+    ! 2020-08-31: bug correction
 
   !=======================================================================
     if (nsboost <= 1) return
@@ -882,11 +884,13 @@ contains
     do i=1, ngr ! loop on low-res rings
        i0 = ringphi(1,i) * nsboost - nsboost - irmin
        do k=-1,1,2 ! West and East side of disc
+          nshrink = 0
           kk = (k+5)/2 ! 2 or 3
 222       continue
           iphi = ringphi(kk, i)
-          if (ringphi(2,i) <= ringphi(3,i) .and. iphi >= 0) then
-             call find_pixel_bounds(nside, nsboost, ringphi(1,i), iphi, phiw, phie)
+          !if (ringphi(2,i) <= ringphi(3,i) .and. iphi >= 0) then
+          if (iphi >= 0) then ! corrected 2020-08-31
+             call find_pixel_bounds(nside, nsboost, ringphi(1,i), iphi, phiw, phie, np=np)
              do j=1, 2*nsboost+1
                 if (i0+j >= 1 .and. i0+j <= nrh) then
                    phic = (phie(j)+phiw(j))*0.5_dp ! pixel center
@@ -896,18 +900,25 @@ contains
                    if (dd <= dph) goto 1000 ! pixel touched by disc, move to next one
                 endif
              enddo
-             ringphi(kk, i)= iphi - k ! pixel not in disc, move edge pixel inwards
-             goto 222 ! try next pixel inward
+             !ringphi(kk, i)= iphi - k ! pixel not in disc, move edge pixel inwards
+             ringphi(kk, i) = mod(iphi + np - k, np) ! pixel not in disc, move edge pixel inwards in [0,np-1]
+             nshrink = nshrink + 1
+             if (nshrink <= 2) then
+                goto 222 ! try next pixel inward
+             else
+                ringphi(kk,i) = -999
+             endif
 1000         continue
           endif
-       enddo ! loop on side
+       enddo ! loop on k East-West sides
     enddo ! loop on low-res rings
 
     ! remove empty rings
     ngr_out = 0
     do i=1,ngr
-       diff = ringphi(3,i) - ringphi(2,i)
-       if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0 .and. diff /= -2 .and. diff /= -1) then
+!        diff = ringphi(3,i) - ringphi(2,i)
+!        if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0 .and. diff /= -2 .and. diff /= -1) then
+       if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0) then ! 2020-08-31
           ngr_out = ngr_out + 1
           ringphi(1:3, ngr_out) = ringphi(1:3, i)
        endif
@@ -921,10 +932,11 @@ contains
     return
   end subroutine check_edge_pixels
   !=======================================================================
-  subroutine find_pixel_bounds (nside, nsboost, iring, iphi, phiw, phie)
+  subroutine find_pixel_bounds (nside, nsboost, iring, iphi, phiw, phie, np)
     !=======================================================================
     integer(i4b),               intent(in)  :: nside, nsboost, iring, iphi
     real(dp),     dimension(1:2*nsboost+1), intent(out) :: phiw, phie
+    integer(i4b), optional,                 intent(out) :: np
     
     real(dp),     dimension(1:2*nsboost+1) :: f, f1, phiw_t, phie_t
     real(dp) :: c0, quad, phie1, phie2, phiw1, phiw2, cv
@@ -936,6 +948,7 @@ contains
     !f = ((/ (i,i=0,2*nsboost) /) - nsboost) / nsboost
     f = ((/ (i*1.d0,i=0,2*nsboost) /) - nsboost*1.d0) / nsboost
 
+    if (present(np)) np = npr
     nq = npr/4 ! number of pixels on current ring in [0,Pi/2] (quadrant)
     transition = (iring == nside .or. iring == nside*3)
 
diff --git a/src/f90/mod/pixel_routines.F90 b/src/f90/mod/pixel_routines.F90
index 26a6b4b..b5d7fac 100644
--- a/src/f90/mod/pixel_routines.F90
+++ b/src/f90/mod/pixel_routines.F90
@@ -1920,6 +1920,7 @@
 ! New algorithm for Inclusive case: test boundary of edge pixels on each ring
 !    2013-04-02: bug correction in query_disc in inclusive mode
 !    2014-07-07: bug correction in discedge2fulldisc
+! 2020-08-31: bug correction in check_edge_pixels (reported by R. de Belsunce)
     !=======================================================================
 #ifdef DOI8B
   subroutine query_disc_8( nside, vector0, radius, listpix, nlist, nest, inclusive)
@@ -2016,8 +2017,9 @@
     call pixels_on_edge(nside, irmin, irmax, phi0, dphitab, ringphi, ngr)
     if (do_inclusive) then
        ! sample edge pixels at larger Nside
-       nsboost = 16
-       nsideh = min(NS_MAX8, nside * int(nsboost,i8b))
+       nsboost = 2**ceiling( log(1.1_DP/(radius*nside))/log(2.0_DP) ) ! make pixel smaller than disc (2020-08-31)
+       nsboost = max(nsboost,16)
+       nsideh  = min(NS_MAX8, nside * int(nsboost,i8b))
        radiush = fudge_query_radius(nsideh, radius, quadratic=.true.)
 
        irmin = ring_num(nsideh, zmax, shift=+1) ! shifted South
@@ -2028,6 +2030,7 @@
           zlist(iz-irmin+1) = ring2z(nsideh, iz)
        enddo
        call discphirange_at_z(vector0, radiush, zlist, nrh, dphilist, phi0)
+       phi0 = mod(phi0 + TWOPI, TWOPI) ! [-Pi,Pi] -> [0,2Pi]
        call check_edge_pixels(nside, nsboost, irmin, irmax, phi0, dphilist, ringphi, ngr)
        deallocate(zlist, dphilist)
     endif
diff --git a/src/f90/mod/sharp_healpix_f_inc.c b/src/f90/mod/sharp_healpix_f_inc.c
index ff7b876..901efb7 100644
--- a/src/f90/mod/sharp_healpix_f_inc.c
+++ b/src/f90/mod/sharp_healpix_f_inc.c
@@ -44,6 +44,7 @@ void X(map2alm_spin) (int nside, int lmax, int mmax, int spin,
   void *mapptr[2], *almptr[2];
 
   CHECK_STACK_ALIGN(8);
+  UTIL_ASSERT(spin>0, "spin must not be zero");
   sharp_make_healpix_geom_info_2 (nside, wgt, zbounds[0], zbounds[1], &ginfo);
   sharp_make_rectangular_alm_info (lmax,mmax,2,&ainfo);
 
@@ -109,6 +110,7 @@ void X(alm2map_spin) (int nside, int lmax, int mmax, int spin,
   void *mapptr[2], *almptr[2];
 
   CHECK_STACK_ALIGN(8);
+  UTIL_ASSERT(spin>0, "spin must not be zero");
   // sharp_make_healpix_geom_info (nside, 1, &ginfo);
   int i,nrings=4*nside-1;
   double *wgt=RALLOC(double,nrings);
diff --git a/src/f90/mod/udgrade_nr.f90 b/src/f90/mod/udgrade_nr.F90
similarity index 97%
rename from src/f90/mod/udgrade_nr.f90
rename to src/f90/mod/udgrade_nr.F90
index b0d6178..757115d 100644
--- a/src/f90/mod/udgrade_nr.f90
+++ b/src/f90/mod/udgrade_nr.F90
@@ -59,10 +59,10 @@ module udgrade_nr
 contains
 
   ! include SP implementation of routines
-  include 'udgrade_s_inc.f90'
+#include "udgrade_s_inc.f90"
   
   ! include DP implementation of routines
-  include 'udgrade_d_inc.f90'
+#include "udgrade_d_inc.f90"
 
 
 end module udgrade_nr