subroutine lalo_to_tile(tf_lalo,mask_lalo,dlats_lalo,dlons_lalo,jdim_lalo,idim_lalo, & tf_tile,mask_tile,xlats_tile,xlons_tile,jdim_tile,idim_tile, & sfcflag,dsearch,miss_fill,bmiss) !-------------------------------------------------------------------------------- ! abstract: Interpolate lon/lat grid to fv3 native grid (tf_lalo => tf_tile) !-------------------------------------------------------------------------------- ! Input ! ! tf_lalo : (idim_lalo,idim_lalo) tf at lat/lon regular grid ! mask_lalo : (idim_lalo,idim_lalo) mask of tf at lat/lon regular grid ! dlats_lalo : (jdim_lalo) latitudes along y direction of lat/lon regular grid points ! dlons_lalo : (idim_lalo) longitudes along x direction of lat/lon regular grid points ! jdim_lalo : number of y dimension of tf_lalo ! idim_lalo : number of x dimension of tf_lalo ! mask_tile : (jdim_tile*idim_tile) mask of tf at cubed sphere grid ! xlats_tile : (jdim_tile*idim_tile) latitudes of all tile grid points ! xlons_tile : (jdim_tile*idim_tile) longitudes of all tile grid points ! jdim_tile : number of y dimension of tf_tile ! idim_tile : number of x dimension of tf_tile ! sfcflag : surface flag (mask type) of the target ! dsearch : maximum search radius in KM ! miss_fill : modes to fill the grids (when failed in search step) ! bmiss : missing value or default value ! ! Output ! ! tf_tile : (jdim_tile*idim_tile) tf at cubed sphere grid implicit none ! input/output real, dimension(idim_lalo,jdim_lalo), intent(in) :: tf_lalo integer, dimension(idim_lalo,jdim_lalo), intent(in) :: mask_lalo real, dimension(jdim_lalo), intent(in) :: dlats_lalo real, dimension(idim_lalo), intent(in) :: dlons_lalo integer, dimension(jdim_tile*idim_tile), intent(in) :: mask_tile real, dimension(jdim_tile*idim_tile), intent(in) :: xlats_tile real, dimension(jdim_tile*idim_tile), intent(in) :: xlons_tile integer, intent(in) :: jdim_lalo,idim_lalo,jdim_tile,idim_tile, & sfcflag,miss_fill real, intent(in) :: dsearch,bmiss real, dimension(jdim_tile*idim_tile), intent(out) :: tf_tile ! Local real, parameter :: deg2rad=3.1415926/180.0 real, dimension(jdim_lalo) :: xlats_lalo real, dimension(idim_lalo) :: xlons_lalo real :: tf,wsum,res_km integer :: itile,jtile integer :: ii,jj,ij,iii,jjj,nintp,mfill,nfill,nsearch,max_search integer :: ilalo,jlalo,ilalop1,jlalop1 integer :: istart,iend,jstart,jend,krad integer, allocatable, dimension(:,:) :: id1,id2,jdc real, allocatable, dimension(:,:,:) :: agrid,s2c print* print*,'interpolate from lat/lon grids to any one grid with known lat/lon' xlats_lalo = dlats_lalo*deg2rad xlons_lalo = dlons_lalo*deg2rad allocate(agrid(idim_tile,jdim_tile,2)) agrid(:,:,1) = reshape (xlons_tile, (/idim_tile,jdim_tile/) ) agrid(:,:,2) = reshape (xlats_tile, (/idim_tile,jdim_tile/) ) agrid = agrid*deg2rad allocate(id1(idim_tile,jdim_tile)) allocate(id2(idim_tile,jdim_tile)) allocate(jdc(idim_tile,jdim_tile)) allocate(s2c(idim_tile,jdim_tile,4)) !---------------------------------------------------------------------- ! compute bilinear weights for each model point from the nearest ! four lalo points. does not account for mask. that ! happens later. !---------------------------------------------------------------------- call remap_coef( 1, idim_tile, 1, jdim_tile, idim_lalo, jdim_lalo, & xlons_lalo, xlats_lalo, id1, id2, jdc, s2c, agrid ) !---------------------------------------------------------------------- !tf_tile will be output. initialize to bmiss. !---------------------------------------------------------------------- tf_tile = bmiss nintp = 0 nsearch = 0 nfill = 0 mfill = 0 !---------------------------------------------------------------------- ! The maximum distance to search is 500 km. how many gaussian ! grid lengths is that? !---------------------------------------------------------------------- res_km = 360.0/real(idim_lalo)*111.0 max_search = ceiling(dsearch/res_km) write(*,'(a,2F7.1,I8)') 'starting ij_loop,dsearch,res_km,max_search ',dsearch,res_km,max_search ij_loop : do ij = 1, jdim_tile*idim_tile !---------------------------------------------------------------------- ! skip non-water points. !---------------------------------------------------------------------- if ( mask_tile(ij) /= 0 ) then cycle ij_loop endif !---------------------------------------------------------------------- ! these are points that are open water !---------------------------------------------------------------------- jtile = (ij-1)/idim_tile + 1 itile = mod(ij,idim_tile) if (itile==0) itile = idim_tile !---------------------------------------------------------------------- ! see if any of the nearest 4 points mask area open water. ! if so, apply tf using bilinear interpolation. !---------------------------------------------------------------------- ilalo = id1(itile,jtile) ilalop1 = id2(itile,jtile) jlalo = jdc(itile,jtile) jlalop1 = jdc(itile,jtile) + 1 if ( mask_lalo(ilalo,jlalo) == sfcflag .or. & mask_lalo(ilalop1,jlalo) == sfcflag .or. & mask_lalo(ilalop1,jlalop1) == sfcflag .or. & mask_lalo(ilalo,jlalop1) == sfcflag ) then tf = 0.0 wsum = 0.0 if (mask_lalo(ilalo,jlalo) == sfcflag) then tf = tf + (s2c(itile,jtile,1)*tf_lalo(ilalo,jlalo)) wsum = wsum + s2c(itile,jtile,1) endif if (mask_lalo(ilalop1,jlalo) == sfcflag) then tf = tf + (s2c(itile,jtile,2)*tf_lalo(ilalop1,jlalo)) wsum = wsum + s2c(itile,jtile,2) endif if (mask_lalo(ilalop1,jlalop1) == sfcflag) then tf = tf + (s2c(itile,jtile,3)*tf_lalo(ilalop1,jlalop1)) wsum = wsum + s2c(itile,jtile,3) endif if (mask_lalo(ilalo,jlalop1) == sfcflag) then tf = tf + (s2c(itile,jtile,4)*tf_lalo(ilalo,jlalop1)) wsum = wsum + s2c(itile,jtile,4) endif if ( wsum > 0.0 ) then nintp = nintp + 1 tf_tile(ij) = tf/wsum else ! write(*,*) 'Warning: No water points nearby, do search/fill' endif else ! handle the case when no nearby open water grids !---------------------------------------------------------------------- ! no nearby gsi/gaussian open water points. perform a spiral search to ! find nearest non-land point on gsi/gaussian grid in distance of ! max_search grid number. !---------------------------------------------------------------------- do krad = 1, max_search istart = ilalo - krad iend = ilalo + krad jstart = jlalo - krad jend = jlalo + krad do jj = jstart, jend do ii = istart, iend if ( (jj == jstart) .or. (jj == jend) .or. & (ii == istart) .or. (ii == iend)) then if ((jj >= 1) .and. (jj <= jdim_lalo)) then jjj = jj if (ii <= 0) then iii = idim_lalo + ii else if (ii >= (idim_lalo+1)) then iii = ii - idim_lalo else iii = ii end if if (mask_lalo(iii,jjj) == sfcflag) then nsearch = nsearch + 1 tf_tile(ij) = tf_lalo(iii,jjj) cycle ij_loop endif ! lalo mask is open water endif endif enddo enddo enddo ! krad loop ! print*,'warning !!!!!! search failed at tile point ',itile,jtile ! ! Fill with nearby 4-point average of SST climatology ! ! write(*,*) 'Warning: fill needed' if ( miss_fill > 0 ) then nfill = nfill + 1 else mfill = mfill + 1 endif endif enddo ij_loop write(*,*) 'nintp = ',nintp write(*,*) 'nsearch = ',nsearch write(*,*) 'nfill = ',nfill write(*,*) 'mfill = ',mfill deallocate(id1, id2, jdc, s2c) end subroutine lalo_to_tile subroutine remap_coef( is, ie, js, je,& im, jm, lon, lat, id1, id2, jdc, s2c, agrid ) !---------------------------------------------------------------------- ! this routine was taken from the forecast model - ! ./atmos_cubed_sphere/tools/fv_treat_da_inc.f90. !---------------------------------------------------------------------- implicit none integer, intent(in) :: is, ie, js, je integer, intent(in) :: im, jm real, intent(in) :: lon(im), lat(jm) real, intent(out) :: s2c(is:ie,js:je,4) integer, intent(out), dimension(is:ie,js:je):: id1, id2, jdc real, intent(in) :: agrid(is:ie,js:je,2) ! local: real :: rdlon(im) real :: rdlat(jm) real :: a1, b1 real, parameter :: pi = 3.1415926 integer i,j, i1, i2, jc, i0, j0 do i=1,im-1 rdlon(i) = 1. / (lon(i+1) - lon(i)) enddo rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im)) do j=1,jm-1 rdlat(j) = 1. / (lat(j+1) - lat(j)) enddo ! * Interpolate to cubed sphere cell center do 5000 j=js,je do i=is,ie if ( agrid(i,j,1)>lon(im) ) then i1 = im; i2 = 1 a1 = (agrid(i,j,1)-lon(im)) * rdlon(im) elseif ( agrid(i,j,1)=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then i1 = i0; i2 = i0+1 a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0) go to 111 endif enddo endif 111 continue if ( agrid(i,j,2)lat(jm) ) then jc = jm-1 b1 = 1. else do j0=1,jm-1 if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then jc = j0 b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc) go to 222 endif enddo endif 222 continue if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then write(*,*) 'gid=', i,j,a1, b1 endif s2c(i,j,1) = (1.-a1) * (1.-b1) s2c(i,j,2) = a1 * (1.-b1) s2c(i,j,3) = a1 * b1 s2c(i,j,4) = (1.-a1) * b1 id1(i,j) = i1 id2(i,j) = i2 jdc(i,j) = jc enddo !i-loop 5000 continue ! j-loop end subroutine remap_coef