program modif
 
USE netcdf
 
IMPLICIT NONE
 
INTEGER :: fidA, status, dimID_x, dimID_y, dimID_time, mx, my, mtime, y_ID, x_ID, smb_ID, sector_ID, time_ID, fidM, &
&          rs, fac, i, j, inew, jnew, mxnew, mynew, iMARclo, jMARclo, iMARinf, iMARsup, jMARinf, jMARsup, nn, si, sj

INTEGER*2,ALLOCATABLE,DIMENSION(:,:) :: msk
 
CHARACTER(LEN=100) :: file_in, file_out
 
REAL*4,ALLOCATABLE,DIMENSION(:) :: y, x, ynew, xnew
 
REAL*4,ALLOCATABLE,DIMENSION(:,:) :: smb, smbnew, xst, yst

REAL*4 :: MAR_lon0, MAR_lat0, MAR_x_stereo_center, MAR_y_stereo_center, theta, deg2rad, xMAReq, yMAReq
 
LOGICAL :: ok

file_in  = 'MARv3.9.3_Amundsen_smb_2088-2117_mean.nc'
file_out = 'MARv3.9.3_Amundsen_smb_2088-2117_mean_polarstereo.nc'
 
MAR_lon0 = -112.0 ! MAR domain center longitude in MARgrd.ctr
MAR_lat0 =  -75.0 ! MAR domain center latitude in MARgrd.ctr

MAR_x_stereo_center = -1519453.35968732 ! from python, in m:  [x,y]=st.lonlat_to_xy([-112],[-75])     
MAR_y_stereo_center =  -613899.00624702 ! from python, in m:  [x,y]=st.lonlat_to_xy([-112],[-75])     

theta = -20.  ! x-Direction (2D runs only ; 0=North, clockwise) in MARgrd.ctr

deg2rad = acos(-1.e0) / 180.e0
write(*,*) 'pi = ', acos(-1.e0)

mxnew = 138
mynew = 140

fac=5

!---------------------------------------
! Read netcdf input file :
 
write(*,*) 'Reading ', TRIM(file_in)
 
status = NF90_OPEN(TRIM(file_in),0,fidA); call erreur(status,.TRUE.,"read")
 
status = NF90_INQ_DIMID(fidA,"x",dimID_x); call erreur(status,.TRUE.,"inq_dimID_x")
status = NF90_INQ_DIMID(fidA,"y",dimID_y); call erreur(status,.TRUE.,"inq_dimID_y")
 
status = NF90_INQUIRE_DIMENSION(fidA,dimID_x,len=mx); call erreur(status,.TRUE.,"inq_dim_x")
status = NF90_INQUIRE_DIMENSION(fidA,dimID_y,len=my); call erreur(status,.TRUE.,"inq_dim_y")
  
ALLOCATE(  y(my), ynew(mynew*fac)  ) 
ALLOCATE(  x(mx), xnew(mxnew*fac)  ) 
ALLOCATE(  smb(mx,my), msk(mx,my), smbnew(mxnew*fac,mynew*fac), xst(mx,my), yst(mx,my)  ) 
 
status = NF90_INQ_VARID(fidA,"y",y_ID); call erreur(status,.TRUE.,"inq_y_ID")
status = NF90_INQ_VARID(fidA,"x",x_ID); call erreur(status,.TRUE.,"inq_x_ID")
status = NF90_INQ_VARID(fidA,"smb",smb_ID); call erreur(status,.TRUE.,"inq_smb_ID")
 
status = NF90_GET_VAR(fidA,y_ID,y); call erreur(status,.TRUE.,"getvar_y")
status = NF90_GET_VAR(fidA,x_ID,x); call erreur(status,.TRUE.,"getvar_x")
status = NF90_GET_VAR(fidA,smb_ID,smb); call erreur(status,.TRUE.,"getvar_smb")
 
status = NF90_CLOSE(fidA); call erreur(status,.TRUE.,"close_file")

where( abs(smb(:,:)) .lt. 1.e9 )
  msk(:,:)=1
elsewhere
  msk(:,:)=0
  smb(:,:) = 1.e9 ! to avoid NaN when multiplying by 0
endwhere
 
!---------------------------------------
! MAR coord in standard stereo:

x(:) = x(:) * 1.e3 ! meters
y(:) = y(:) * 1.e3 ! meters

theta = theta * deg2rad

do i=1,mx
do j=1,my
  xst(i,j) = MAR_x_stereo_center + x(i)*sin(theta) - y(j)*cos(theta)
  yst(i,j) = MAR_y_stereo_center + x(i)*cos(theta) + y(j)*sin(theta)
enddo
enddo  

write(*,*) MINVAL(xst), MAXVAL(xst)
write(*,*) MINVAL(yst), MAXVAL(yst)

! new grid

do inew=1,mxnew*fac
 xnew(inew) = -2.e6 + (inew-1-fac/2)*8000./fac ! meters
enddo

do jnew=1,mynew*fac
 ynew(jnew) = -1.e6 + (jnew-1-fac/2)*8000./fac ! meters
enddo

! interpolation:
do inew=1,mxnew*fac
do jnew=1,mynew*fac

  xMAReq = ( ( xnew(inew) - MAR_x_stereo_center ) * sin(theta) + ( ynew(jnew) - MAR_y_stereo_center ) * cos(theta) )
  yMAReq = ( ( MAR_x_stereo_center - xnew(inew) ) * cos(theta) + ( ynew(jnew) - MAR_y_stereo_center ) * sin(theta) )

  iMARclo = MINLOC( ABS( x - xMAReq ),1 )
  jMARclo = MINLOC( ABS( y - yMAReq ),1 )

  if ( ABS(x(MIN(iMARclo+1,mx))-xnew(inew)) .gt. ABS(x(MAX(iMARclo-1,1))-xnew(inew)) ) then
    iMARinf=MAX(iMARclo-1,1)
    iMARsup=iMARclo
  else
    iMARinf=iMARclo
    iMARsup=MIN(iMARclo+1,mx)
  endif
  if ( ABS(y(MIN(jMARclo+1,my))-ynew(jnew)) .gt. ABS(y(MAX(jMARclo-1,1))-ynew(jnew)) ) then
    jMARinf=MAX(jMARclo-1,1)
    jMARsup=jMARclo
  else
    jMARinf=jMARclo
    jMARsup=MIN(jMARclo+1,my)
  endif

  if ( msk(iMARinf,jMARinf) + msk(iMARinf,jMARsup) + msk(iMARsup,jMARinf) + msk(iMARsup,jMARsup) .ge. 1 ) then  
    smbnew(inew,jnew) =   (   abs(0.5*(xst(iMARsup,jMARinf)+xst(iMARsup,jMARsup))-xnew(inew)) * abs(0.5*(yst(iMARsup,jMARinf)+yst(iMARinf,jMARinf))-ynew(jnew)) * smb(iMARinf,jMARinf) * msk(iMARinf,jMARinf)  &
    &                       + abs(0.5*(xst(iMARsup,jMARinf)+xst(iMARsup,jMARsup))-xnew(inew)) * abs(0.5*(yst(iMARsup,jMARsup)+yst(iMARinf,jMARsup))-ynew(jnew)) * smb(iMARinf,jMARsup) * msk(iMARinf,jMARsup)  & 
    &                       + abs(0.5*(xst(iMARinf,jMARinf)+xst(iMARinf,jMARsup))-xnew(inew)) * abs(0.5*(yst(iMARsup,jMARinf)+yst(iMARinf,jMARinf))-ynew(jnew)) * smb(iMARsup,jMARinf) * msk(iMARsup,jMARinf)  &
    &                       + abs(0.5*(xst(iMARinf,jMARinf)+xst(iMARinf,jMARsup))-xnew(inew)) * abs(0.5*(yst(iMARsup,jMARsup)+yst(iMARinf,jMARsup))-ynew(jnew)) * smb(iMARsup,jMARsup) * msk(iMARsup,jMARsup) )& 
    &                   * 1.e-3 * 365.25  & ! mm.w.e/day -> m.w.e /yr 
    &                   / (   abs(0.5*(xst(iMARsup,jMARinf)+xst(iMARsup,jMARsup))-xnew(inew)) * abs(0.5*(yst(iMARsup,jMARinf)+yst(iMARinf,jMARinf))-ynew(jnew))                        * msk(iMARinf,jMARinf)  &
    &                       + abs(0.5*(xst(iMARsup,jMARinf)+xst(iMARsup,jMARsup))-xnew(inew)) * abs(0.5*(yst(iMARsup,jMARsup)+yst(iMARinf,jMARsup))-ynew(jnew))                        * msk(iMARinf,jMARsup)  & 
    &                       + abs(0.5*(xst(iMARinf,jMARinf)+xst(iMARinf,jMARsup))-xnew(inew)) * abs(0.5*(yst(iMARsup,jMARinf)+yst(iMARinf,jMARinf))-ynew(jnew))                        * msk(iMARsup,jMARinf)  &
    &                       + abs(0.5*(xst(iMARinf,jMARinf)+xst(iMARinf,jMARsup))-xnew(inew)) * abs(0.5*(yst(iMARsup,jMARsup)+yst(iMARinf,jMARsup))-ynew(jnew))                        * msk(iMARsup,jMARsup) ) 

    !smbnew(inew,jnew) =   (   abs(xst(iMARsup,jMARsup)-xnew(inew)) * abs(yst(iMARsup,jMARsup)-ynew(jnew)) * smb(iMARinf,jMARinf) * msk(iMARinf,jMARinf)   &
    !&                       + abs(xst(iMARsup,jMARinf)-xnew(inew)) * abs(ynew(jnew)-yst(iMARsup,jMARinf)) * smb(iMARinf,jMARsup) * msk(iMARinf,jMARsup)   &
    !&                       + abs(xnew(inew)-xst(iMARinf,jMARsup)) * abs(yst(iMARinf,jMARsup)-ynew(jnew)) * smb(iMARsup,jMARinf) * msk(iMARsup,jMARinf)   &
    !&                       + abs(xnew(inew)-xst(iMARinf,jMARinf)) * abs(ynew(jnew)-yst(iMARinf,jMARinf)) * smb(iMARsup,jMARsup) * msk(iMARsup,jMARsup) ) &
    !&                   * 1.e-3 * 365.25  & ! mm.w.e/day -> m.w.e /yr 
    !&                   / (   abs(xst(iMARsup,jMARsup)-xnew(inew)) * abs(yst(iMARsup,jMARsup)-ynew(jnew))                        * msk(iMARinf,jMARinf)   &
    !&                       + abs(xst(iMARsup,jMARinf)-xnew(inew)) * abs(ynew(jnew)-yst(iMARsup,jMARinf))                        * msk(iMARinf,jMARsup)   &
    !&                       + abs(xnew(inew)-xst(iMARinf,jMARsup)) * abs(yst(iMARinf,jMARsup)-ynew(jnew))                        * msk(iMARsup,jMARinf)   &
    !&                       + abs(xnew(inew)-xst(iMARinf,jMARinf)) * abs(ynew(jnew)-yst(iMARinf,jMARinf))                        * msk(iMARsup,jMARsup) )

  else
    smbnew(inew,jnew) = NF90_FILL_FLOAT
  endif

  !write(*,*) '##### ', inew, jnew
  !write(*,*) '      ', iMARinf, iMARsup, jMARinf, jMARsup
  !write(*,*) '      ', smb(iMARinf,jMARinf), smb(iMARinf,jMARsup), smb(iMARsup,jMARinf), smb(iMARsup,jMARsup)
  !write(*,*) '      ', smbnew(inew,jnew)

  if ( inew .eq. 9 .and. jnew .eq. 497 ) then
      write(*,*) inew,jnew,smbnew(inew,jnew)
      write(*,*) iMARinf, iMARsup, jMARinf, jMARsup
      write(*,*) abs(0.5*(xst(iMARsup,jMARinf)+xst(iMARsup,jMARsup))-xnew(inew)) , abs(0.5*(yst(iMARsup,jMARinf)+yst(iMARinf,jMARinf))-ynew(jnew)) , smb(iMARinf,jMARinf) , msk(iMARinf,jMARinf)
      write(*,*) abs(0.5*(xst(iMARsup,jMARinf)+xst(iMARsup,jMARsup))-xnew(inew)) , abs(0.5*(yst(iMARsup,jMARsup)+yst(iMARinf,jMARsup))-ynew(jnew)) , smb(iMARinf,jMARsup) , msk(iMARinf,jMARsup)
      write(*,*) abs(0.5*(xst(iMARinf,jMARinf)+xst(iMARinf,jMARsup))-xnew(inew)) , abs(0.5*(yst(iMARsup,jMARinf)+yst(iMARinf,jMARinf))-ynew(jnew)) , smb(iMARsup,jMARinf) , msk(iMARsup,jMARinf)
      write(*,*) abs(0.5*(xst(iMARinf,jMARinf)+xst(iMARinf,jMARsup))-xnew(inew)) , abs(0.5*(yst(iMARsup,jMARsup)+yst(iMARinf,jMARsup))-ynew(jnew)) , smb(iMARsup,jMARsup) , msk(iMARsup,jMARsup)
  endif


enddo
enddo

! re-average over coarser grid
DEALLOCATE(x, y, smb)
ALLOCATE( x(mxnew), y(mynew), smb(mxnew,mynew) )

do inew=1,mxnew
 x(inew) = SUM( xnew(fac*(inew-1)+1:fac*(inew-1)+fac) ) / fac
enddo

do jnew=1,mynew
 y(jnew) = SUM( ynew(fac*(jnew-1)+1:fac*(jnew-1)+fac) ) / fac
enddo

smb(:,:) = 0.e0
do inew=1,mxnew
do jnew=1,mynew
 nn=0
 do si=1,fac
 do sj=1,fac
   if ( smbnew(MAX(MIN(fac*(inew-1)+si,mxnew*fac),1),MAX(MIN(fac*(jnew-1)+sj,mynew*fac),1)) .ne. NF90_FILL_FLOAT ) then
     smb(inew,jnew) = smb(inew,jnew) + smbnew(MAX(MIN(fac*(inew-1)+si,mxnew*fac),1),MAX(MIN(fac*(jnew-1)+sj,mynew*fac),1))
     nn=nn+1
   endif
 enddo
 enddo
 if ( nn .gt. 0 ) then
     !write(*,*) smb(inew,jnew), nn
     smb(inew,jnew) = smb(inew,jnew) / nn
 else
     smb(inew,jnew) = NF90_FILL_FLOAT
 endif
enddo
enddo

! extrapolate to closest locations with missing values (in case the ice-sheet extent i s a bit larger)
DEALLOCATE( smbnew )
ALLOCATE( smbnew(mxnew,mynew) )

do rs=1,100
  smbnew(:,:) = smb(:,:)
  do inew=1,mxnew
  do jnew=1,mynew
    if ( smb(inew,jnew) .eq. NF90_FILL_FLOAT ) then
      !write(*,*) '%% ', rs, inew,jnew
      ok=.false.
      do si=-1,1
        do sj=-1,1
          !write(*,*) si, sj, smb(MAX(MIN(inew+si,mxnew),1),MAX(MIN(jnew+sj,mynew),1)), ok 
          if ( smb(MAX(MIN(inew+si,mxnew),1),MAX(MIN(jnew+sj,mynew),1)) .ne. NF90_FILL_FLOAT ) then
            !write(*,*) '   ', smbnew(inew,jnew), smb(MAX(MIN(inew+si,mxnew),1),MAX(MIN(jnew+sj,mynew),1))
            smbnew(inew,jnew) = smb(MAX(MIN(inew+si,mxnew),1),MAX(MIN(jnew+sj,mynew),1))
            ok=.true.
            !write(*,*) '   ', MAX(MIN(inew+si,mxnew),1),MAX(MIN(jnew+sj,mynew),1)
          endif
          if (ok) exit
        enddo
        if (ok) exit
      enddo
    endif
  enddo
  enddo
  smb(:,:) = smbnew(:,:)
enddo

!---------------------------------------
! Writing new netcdf file :
 
write(*,*) 'Creating ', TRIM(file_out)
 
status = NF90_CREATE(TRIM(file_out),NF90_NOCLOBBER,fidM); call erreur(status,.TRUE.,'create')
 
status = NF90_DEF_DIM(fidM,"x",mxnew,dimID_x); call erreur(status,.TRUE.,"def_dimID_x")
status = NF90_DEF_DIM(fidM,"y",mynew,dimID_y); call erreur(status,.TRUE.,"def_dimID_y")
  
status = NF90_DEF_VAR(fidM,"y",NF90_FLOAT,(/dimID_y/),y_ID); call erreur(status,.TRUE.,"def_var_y_ID")
status = NF90_DEF_VAR(fidM,"x",NF90_FLOAT,(/dimID_x/),x_ID); call erreur(status,.TRUE.,"def_var_x_ID")
status = NF90_DEF_VAR(fidM,"smb",NF90_FLOAT,(/dimID_x,dimID_y/),smb_ID); call erreur(status,.TRUE.,"def_var_smb_ID")
 
status = NF90_PUT_ATT(fidM,y_ID,"standard_name","y"); call erreur(status,.TRUE.,"put_att_y_ID")
status = NF90_PUT_ATT(fidM,y_ID,"long_name","y"); call erreur(status,.TRUE.,"put_att_y_ID")
status = NF90_PUT_ATT(fidM,y_ID,"units","m"); call erreur(status,.TRUE.,"put_att_y_ID")
status = NF90_PUT_ATT(fidM,x_ID,"standard_name","x"); call erreur(status,.TRUE.,"put_att_x_ID")
status = NF90_PUT_ATT(fidM,x_ID,"long_name","x"); call erreur(status,.TRUE.,"put_att_x_ID")
status = NF90_PUT_ATT(fidM,x_ID,"units","m"); call erreur(status,.TRUE.,"put_att_x_ID")
status = NF90_PUT_ATT(fidM,smb_ID,"standard_name","Surface_Mass_Balance"); call erreur(status,.TRUE.,"put_att_smb_ID")
status = NF90_PUT_ATT(fidM,smb_ID,"long_name","Surface Mass Balance"); call erreur(status,.TRUE.,"put_att_smb_ID")
status = NF90_PUT_ATT(fidM,smb_ID,"units","m.w.e/yr"); call erreur(status,.TRUE.,"put_att_smb_ID")
 
status = NF90_PUT_ATT(fidM,NF90_GLOBAL,"history","Created using interpolate_to_std_polar_stereographic.f90"); call erreur(status,.TRUE.,"put_att_GLOBAL_ID")
 
status = NF90_ENDDEF(fidM); call erreur(status,.TRUE.,"fin_definition") 
 
status = NF90_PUT_VAR(fidM,y_ID,y); call erreur(status,.TRUE.,"var_y_ID")
status = NF90_PUT_VAR(fidM,x_ID,x); call erreur(status,.TRUE.,"var_x_ID")
status = NF90_PUT_VAR(fidM,smb_ID,smb); call erreur(status,.TRUE.,"var_smb_ID")
 
! Close file
status = NF90_CLOSE(fidM); call erreur(status,.TRUE.,"final")

end program modif



SUBROUTINE erreur(iret, lstop, chaine)
  ! pour les messages d'erreur
  USE netcdf
  INTEGER, INTENT(in)                     :: iret
  LOGICAL, INTENT(in)                     :: lstop
  CHARACTER(LEN=*), INTENT(in)            :: chaine
  !
  CHARACTER(LEN=80)                       :: message
  !
  IF ( iret .NE. 0 ) THEN
    WRITE(*,*) 'ROUTINE: ', TRIM(chaine)
    WRITE(*,*) 'ERREUR: ', iret
    message=NF90_STRERROR(iret)
    WRITE(*,*) 'CA VEUT DIRE:',TRIM(message)
    IF ( lstop ) STOP
  ENDIF
  !
END SUBROUTINE erreur
