Problem: Tiling in RSM-> ROMS creates too big steps in values at one neighboring points and zero variations in other neighboring points.
This generates a huge noice in curl and divergence of wind stresses...
I may have to just use gridinterp.f90/..

RSM -> ROMS: tiling.f
ROMS -> RSM: tile_averge.f

October 1, 2007
Hyodae Seo

Instead of using time-consuming gridinterp.f90 and fillnan.f90, why not just simply tiling or averge over the different grid box during the course of interpolation?

RSM grid size=100km * 100k
20
->
ROMS grid size=25km * 25k
20
20
20
20
20
20
20
20
20
20
20
20
20
20
20
20

In order to use these, you have to create the grid for RSM and ROMS in such a manner that,
1 RSM grid box (usually coraser than ROMS, say here 100km resolution ) covers 16 grid boxes of ROMS.
To do that, you have to modify the make_grid.m file for ROMS.
Create RSM grid first and get the grid information from ctl file and modify the portaion of the matlat code that computes lonr and latr.

% use exact 1/4 grids of RSM
clear a;clear b
a.lonr=[111.038:1.038:251.1680];
        a.dlon=1.038;
a.latr=[ 18.985  19.964  20.937  21.903  22.864  23.817  24.763  25.703  26.635  27.559 ...
  28.476  29.385  30.285  31.178  32.062  32.938  33.805  34.664  35.513  36.354 ...
  37.186  38.009  38.822  39.627  40.422  41.208  41.984  42.752  43.509  44.258 ...
  44.997  45.726  46.447  47.157  47.859  48.551  49.234  49.907  50.571  51.226 ...
  51.872  52.508  53.136  53.754  54.364  54.964  55.556  56.139  56.713  57.279 ...
  57.836  58.385  58.925  59.457  59.980  60.496  61.003  61.502  61.994  62.478 ...
  62.954  63.422  63.883  64.336  64.782  65.221];
 
        for i=1:length(a.latr)-1;
        a.dlat1(i)=(a.latr(i+1)-a.latr(i));
        end;
        a.dlat=a.dlat1/4;
b.latr=zeros(1,65*3+66);
b.latr(1:4:end)=a.latr;
ii=0;
for i=1:65;disp(i);
ii=ii+1;
index=i*4-2:i*4;disp(index);
b.latr(index)=[a.latr(i)+a.dlat(ii)*1 a.latr(ii)+a.dlat(ii)*2  a.latr(ii)+a.dlat(ii)*3 ];
end;
 
b.lonr=[a.lonr(1):a.dlon/4:a.lonr(end)];
 
lonr=b.lonr;
latr=b.latr;
 
This code just uniformly divide the 1 grid box of RSM into 4 grids in x and y direction, yielding 16 grid box for ROMS.

RSM to ROMS: tiling.f does just produce a tile (4 by 4) of RSM output that covers 16 ROMS grids.

ROMS grid size=25km * 25k
land
land
20
22
land
land
20
23
land
25
23
20
land
27
26
20
-->

RSM grid size=100km * 100k
(20+22+20+23+25+23+20+27+
26+20 )/10=22.6

tiling.f
      real, dimension(:,:), allocatable :: datain
      real, dimension(:,:), allocatable :: dataout
 
! read input data
          open(11,file='fort.11',form='formatted')
          read(11,*) nx1, ny1
          allocate(datain(nx1,ny1))
 
! read input grid info
          open(12,file='fort.12',form='formatted')
          read(12,*) nx2, ny2
          allocate(dataout(nx2,ny2))
 
! output grid after tiling
         open (unit=13,file='fort.13', form='formatted')
         read (13,*) datain
 
! tiling
          do i1=1,nx1-1
           do j1=1,ny1-1
                  do ii=(i1-1)*4+1,(i1-1)*4+4
                  do jj=(j1-1)*4+1,(j1-1)*4+4
                    dataout(ii,jj)=datain(i1,j1)
                  enddo
                  enddo
             enddo
            enddo
 
! at the northern edge
                do i1=1,nx1-1
                  do ii=(i1-1)*4+1,(i1-1)*4+4
                    dataout(ii,ny2)=datain(i1,ny1)
                  enddo
               enddo
 
! at the eastern edge
           do j1=1,ny1-1
                  do jj=(j1-1)*4+1,(j1-1)*4+4
                    dataout(nx2,jj)=datain(nx1,j1)
                  enddo
             enddo
 
! edge point of at nx2 and ny2
 
            dataout(nx2,ny2)=datain(nx1,ny1)
 
     open (unit=51,file='fort.51', form='formatted')
          write(51,*) dataout
      stop
      end
 
 

ROMS to RSM: tile_averge.f averages SST over the 16 grids box to produce 1 SST for the corresponding RSM grid.
Since land-sea mask is not identical, I clumsily did in fortran something like nanmean in matlab.
Also note that there ARE STILL INCONSISTENT GRID POINTS WHERE SST ARE 0 (or 273.15).
After checking the SST grib file, I manually fill those bad points with the manually selected neighboring point..
See below the cods..
For the diffferent domain, this needs to to modified.. Therefore, tiel_average needs to be created for each domain.
For NP4 case, ile_averge_NP4.f is used..

tile_average_NP4.f
      parameter(rland=99999.)
      real, dimension(:,:), allocatable :: datain
      real, dimension(:,:), allocatable :: dataout
      real, dimension(:,:), allocatable :: rmask
 
! read input data
          open(11,file='fort.11',form='formatted')
          read(11,*) nx1, ny1, nx1ny1
          allocate(datain(nx1,ny1))
          allocate(rmask(nx1,ny1))
 
! read input grid info
          open(12,file='fort.12',form='formatted')
          read(12,*) nx2, ny2
          allocate(dataout(nx2,ny2))
 
! read input grid info
          open(13,file='fort.13',form='formatted')
          read(13,*) rmask
 
! output grid after tiling
         open (unit=14,file='fort.14', form='formatted')
         read (14,*) datain
 
! insert nan value at land points.
         do i=1,nx1
            do j=1,ny1
                   if(rmask(i,j) .eq. 0 ) then
                      datain(i,j)=rland
                   endif
                 enddo
                 enddo
 
! tiling
          do i1=1,nx2-1
           do j1=1,ny2-1
                  summ=0.
                  number_land=0
                  number_sea=0
 
                  do ii=(i1-1)*4+1,(i1-1)*4+4
                  do jj=(j1-1)*4+1,(j1-1)*4+4
 
                     if(datain(ii,jj) .eq. rland) then
                      number_land=number_land+1
                      summ=summ+0
                     else
                      number_sea=number_sea+1
                     summ=summ+datain(ii,jj)
                      endif
 
                  enddo
                  enddo
 
! quick and dirty way of doing nanmean
                    number_landsea=number_land+number_sea
                    if (number_landsea .ne. 16) then
                        print *, "number_land ~= number_sea"
                        STOP
                    endif
 
                     if (number_sea .eq. 0) then
! all points are land, put sst=0;
                        dataout(i1,j1)=0
                     else
                        dataout(i1,j1)=summ/number_sea
                     endif
             enddo
            enddo
 
! nx1,ny1,datain=large
! nxy,ny2,dataout=small
 
 nx1,ny1,datain=large
! nxy,ny2,dataout=small
 
! nothern and eastern boundary
! just read the neighboring points
! nothern:
!             dataout(i2,ny2)=dataout(i2,ny2-1)
! eastern
!             dataout(nx2,j2)=dataout(nx2-1,j2)
 
 
! at the northern edge
          do i2=1,nx2
            dataout(i2,ny2)=dataout(i2,ny2-1)
          enddo
! at the eastern edge
          do j2=1,ny2
            dataout(nx2,j2)=dataout(nx2-1,j2)
          enddo
 
 
 
! find bad SST values and fill them with some other "reasonable and meaningful" values...
! you first run the model without this and manually find wrong SSTs (usually near the coast)
! SST are usually overly or somewhat low than the neighboring points becuase it interpolated with the SST
! over the land grid point where SSTs are zero.
 
 !               near Russia
                  dataout(45,56)= dataout(44,56)
!               near west of Hokkaido
                  dataout(28,29)=dataout(29,29)
!               near Alaska
                  dataout(78,47)=dataout(77,47)
!                 lake near SF
                  dataout(113,36)=0.
!                 lake in Canada
                  dataout(120,56)=0.
                  dataout(120,57)=0.
                  dataout(120,58)=0.
                  dataout(120,59)=0.
 
 
                  dataout(5,8)=dataout(5,7)
                  dataout(6,9)=dataout(7,9)
 
 
     open (unit=51,file='fort.51', form='formatted')
          write(51,*) dataout
      stop
      end
 
 
These will replace gridinterp.f90 in Rsm2Roms.sh and Roms2Rsm.sh.

Things to do: make it general (currently only for tiling is 16).