
c        1         2         3         4         5         6         7
c2345>789012345678901234567890123456789012345678901234567890123456789012

cI    The subroutine "half_resolution" produces an output field with
cI    half the resolution of the input field. The half resolution field
cI    may cover the full domain of the original field, or a subdomain.
cI    The algorithm used ensures that all original nodes are assigned
cI    with equal weight, with the exception of those that are next to
cI    undefined (dry) nodes or the edges. The subroutine may be called
cI    twice to produce one quarter resolution output etc. Remember to
cI    reset 'new_map' to '.true.' everytime a new land-sea grid is used,
cI    e.g., when an original half resolution output is processed to a
cI    quarter resolution product.

cIO   Input:
cIO    nx_in,
cIO    ny_in    (i) - dimension of the fine mesh (original) domain,
cIO                   used for 'field_in(,)'
cIO    field_in
cIO            (RR) - the original fine mesh field
cIO    x0,y0    (i) - fine mesh coordinates of the lower left corner
cIO                   of the subdomain to consider
cIO    xN,yN    (i) - fine mesh coordinates of the upper right corner
cIO                   of the subdomain to consider
cIO    new_map  (l) - should be set to '.true.' when 'map_factor(,,)'
cIO                   must be set (or reset), i.e when a new land-sea
cIO                   grid is encounted.
cIO    nx_out,
cIO    ny_out   (i) - dimension of the coarse mesh (output) domain,
cIO                   used for 'field_out(,)', 'map_factors(,,)' and
cIO                   'defined_node(,)'

cIO   Input/ output:
cIO    map_factor(9,,)
cIO           (RRR) - facors that will be used in the calculation of
cIO                   the coarse grid field, set when 'new_map' is
cIO                   '.true.', should not be modified outside this
cIO                   subroutine
cIO    defined_node(,)
cIO            (LL) - specification of coarse mesh land-sea grid, set
cIO                   when 'new_map' is '.true.', should not be
cIO                   modified outside this subroutine

cIO   Output:
cIO    field_out
cIO            (RR) - the coarse mesh output field

cE    Example:
cE     Input:
cE               nx_in = ny_in = 8
cE               x0 = y0 = 1
cE               xN = yN = 8
cE               new_map = .true.
cE               nx_out = ny_out = 4

cE        Example 1:

cE        Example 2:

cFS   Dependencies:
cFS    External: none
cFS    Internal: functions get_value, get_edge_x, get_edge_y

cR    2000-02-15 (AM):
cR     Corrected error in the two final if-blocks of "half_resolution"


      subroutine half_resolution(nx_in,ny_in,field_in,x0,y0,xN,yN,
     >                           new_map,map_factor,defined_node,
     >                           nx_out,ny_out,field_out)

      implicit none

c     Input:
      integer nx_in,ny_in,x0,y0,xN,yN,nx_out,ny_out
      real field_in(nx_in,ny_in)

c     Input/ output:
      logical new_map,defined_node(nx_out,ny_out)
      real map_factor(9,nx_out,ny_out)

c     Output:
      real field_out(nx_out,ny_out)

c     Parameters:
      real limit,undef_out
      parameter (limit= 100000., undef_out= -1.e10)

c     Local:
      integer x,y,x_top,y_top,x_coarse,y_coarse,x_fine,y_fine,
     >        x_last,y_last,i,weight

c     Function:
      real get_value,get_edge_x,get_edge_y


c     Check dimensions:
      if ((xN-x0+1)/2 .ne. nx_out .or. (yN-y0+1)/2 .ne. ny_out)
     >   go to 999

c     Set factors?
      if (new_map) then

c        Lower left coarse mesh corner, in fine mesh grids: (x0+1,y0+1)
c        Upper right coarse mesh corner, in fine mesh grids:
         x_top = xN
         if (mod(xN - x0, 2) .eq. 0) x_top = x_top -1
         y_top = yN
         if (mod(yN - y0, 2) .eq. 0) y_top = y_top -1

         y_coarse = 0
         do y = y0+1,y_top,2
            y_coarse = y_coarse + 1
            x_coarse = 0
            do x = x0+1,x_top,2
               x_coarse = x_coarse + 1

c                                                    7  6  5
c              Relative numbering in fine grid,
c                                                    8  9  4
c              where '9' is a coarse grid node:
c                                                    1  2  3


c                                                    
c              The weight function is defined        1  2  1
c              so that each node
c              has the same total weight (with       2  4  2
c              the exception of coastal nodes):
c                                                    1  2  1

               weight = 0
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x-1, y-1, 1, 1)
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x,   y-1, 2, 2)
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x+1, y-1, 3, 1)
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x+1, y  , 4, 2)
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x+1, y+1, 5, 1)
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x  , y+1, 6, 2)
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x-1, y+1, 7, 1)
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x-1, y  , 8, 2)
               call update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                         map_factor,weight,x_coarse,y_coarse,
     >                         x  , y  , 9, 4)
               if (weight .lt. 8) then
                  defined_node(x_coarse,y_coarse) = .false.
               else
                  do i = 1,9
                     map_factor(i,x_coarse,y_coarse) =
     >                  map_factor(i,x_coarse,y_coarse)/float(weight)
                  end do
                  defined_node(x_coarse,y_coarse) = .true.
               end if

            end do
         end do

      end if

c      write(*,*)(map_factor(x,1,1), x=1,9)

c     Compute 'field_out()':

      x_last = nx_out
      if (x_top .eq. nx_in) x_last = x_last - 1
      y_last = ny_out
      if (y_top .eq. ny_in) y_last = y_last - 1

      y_fine = y0 - 1
      do y = 1,y_last
         x_fine = x0 - 1
         y_fine = y_fine + 2

         do x = 1,x_last
            x_fine = x_fine + 2

            if (defined_node(x,y)) then
               field_out(x,y) = get_value(x,y,nx_in,ny_in,field_in,
     >                                    x_fine,y_fine,
     >                                    nx_out,ny_out,map_factor)
            else
               field_out(x,y) = undef_out
            end if

         end do
      end do

      if (x_top .eq. nx_in) then

         y_fine = y0 - 1
         do y = 1,ny_out
            y_fine = y_fine + 2
            if (defined_node(nx_out,y)) then
               field_out(nx_out,y) = get_edge_x(y,nx_in,ny_in,field_in,
     >                                          y_fine,nx_out,ny_out,
     >                                          map_factor)
            else
               field_out(nx_out,y) = undef_out
            end if
         end do

      end if

      if (y_top .eq. ny_in) then

         x_fine = x0 - 1
         do x = 1,nx_out
               x_fine = x_fine + 2
            if (defined_node(x,ny_out)) then
               field_out(x,ny_out) = get_edge_y(x,nx_in,ny_in,field_in,
     >                                          x_fine,nx_out,ny_out,
     >                                          map_factor)
            else
               field_out(x,ny_out) = undef_out
            end if
         end do

      end if


      go to 1000

 999  write(*,*)
      write(*,*)'Error in the subroutine "half_resolution".'
      write(*,*)' Error indimensions:'
      write(*,*)'  Fine mesh: x from ',x0,' to ',xN,' -> ',(xN-x0+1)/2,
     >          ' coarse mesh grids'
      write(*,*)'  Specified coarse mesh grids/ x: ',nx_out
      write(*,*)'  Fine mesh: y from ',y0,' to ',yN,' -> ',(yN-y0+1)/2,
     >          ' coarse mesh grids'
      write(*,*)'  Specified coarse mesh grids/ y: ',ny_out
      write(*,*)'Terminates.'
      write(*,*)
      stop

 1000 return
      end



      subroutine update_map(limit,nx_in,ny_in,field_in,nx_out,ny_out,
     >                      map_factor,weight,x_coarse,y_coarse,
     >                      x,y,position,factor)

      implicit none

c     Input:
      integer nx_in,ny_in,nx_out,ny_out,x,y,x_coarse,y_coarse,
     >        position,factor
      real limit,field_in(nx_in,ny_in)

c     Input/ output:
      integer weight
      real map_factor(9,nx_out,ny_out)


      if (x .gt. nx_in .or. y .gt. ny_in) then
         map_factor(position,x_coarse,y_coarse) = 0.
      else if (abs(field_in(x,y)) .lt. limit) then
         map_factor(position,x_coarse,y_coarse) = float(factor)
         weight = weight + factor
      else
         map_factor(position,x_coarse,y_coarse) = 0.
      end if

      return
      end



      real function get_value(x,y,nx_in,ny_in,field_in,x_fine,y_fine,
     >                        nx_out,ny_out,map_factor)

      implicit none

c     Input:
      integer x,y,nx_in,ny_in,nx_out,ny_out,x_fine,y_fine
      real field_in(nx_in,ny_in),map_factor(9,nx_out,ny_out)

      get_value = map_factor(1,x,y)*field_in(x_fine-1,y_fine-1) +
     >            map_factor(2,x,y)*field_in(x_fine,  y_fine-1) +
     >            map_factor(3,x,y)*field_in(x_fine+1,y_fine-1) +
     >            map_factor(4,x,y)*field_in(x_fine+1,y_fine)   +
     >            map_factor(5,x,y)*field_in(x_fine+1,y_fine+1) +
     >            map_factor(6,x,y)*field_in(x_fine,  y_fine+1) +
     >            map_factor(7,x,y)*field_in(x_fine-1,y_fine+1) +
     >            map_factor(8,x,y)*field_in(x_fine-1,y_fine)   +
     >            map_factor(9,x,y)*field_in(x_fine,  y_fine)

      end



      real function get_edge_x(y,nx_in,ny_in,field_in,y_fine,
     >                         nx_out,ny_out,map_factor)

      implicit none

c     Input:
      integer y,nx_in,ny_in,nx_out,ny_out,y_fine
      real field_in(nx_in,ny_in),map_factor(9,nx_out,ny_out)

      get_edge_x = map_factor(1,nx_out,y)*field_in(nx_in-1,y_fine-1) +
     >             map_factor(2,nx_out,y)*field_in(nx_in,  y_fine-1) +
     >             map_factor(6,nx_out,y)*
     >              field_in(nx_in,  min(y_fine+1,ny_in)) +
     >             map_factor(7,nx_out,y)*
     >              field_in(nx_in-1,min(y_fine+1,ny_in)) +
     >             map_factor(8,nx_out,y)*field_in(nx_in-1,y_fine)   +
     >             map_factor(9,nx_out,y)*field_in(nx_in,  y_fine)

      end



      real function get_edge_y(x,nx_in,ny_in,field_in,x_fine,
     >                         nx_out,ny_out,map_factor)

      implicit none

c     Input:
      integer x,nx_in,ny_in,nx_out,ny_out,x_fine
      real field_in(nx_in,ny_in),map_factor(9,nx_out,ny_out)

      get_edge_y = map_factor(1,x,ny_out)*field_in(x_fine-1,ny_in-1) +
     >             map_factor(2,x,ny_out)*field_in(x_fine,  ny_in-1) +
     >             map_factor(3,x,ny_out)*
     >              field_in(min(x_fine+1,nx_in),ny_in-1) +
     >             map_factor(4,x,ny_out)*
     >              field_in(min(x_fine+1,nx_in),ny_in)   +
     >             map_factor(8,x,ny_out)*field_in(x_fine-1,ny_in)   +
     >             map_factor(9,x,ny_out)*field_in(x_fine,ny_in)

      end
