LCOV - code coverage report
Current view: | top level - chemistry/utils - horizontal_interpolate.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 96 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 2 | 0.0 %

              Line data    Source code
           1             : module horizontal_interpolate
           2             :
           3             : !
           4             : ! Modules Used:
           5             : !
           6             :   use shr_kind_mod,       only: r8 => shr_kind_r8
           7             :   use shr_const_mod,      only: SHR_CONST_PI
           8             :   use cam_abortutils,     only: endrun
           9             :   use scamMod,            only: single_column
          10             :   use cam_logfile,        only: iulog
          11             :   implicit none
          12             :   private
          13             :   save
          14             :
          15             :   real(r8) :: gw1(1000), gw2(1000)
          16             :
          17             :   public :: xy_interp_init, xy_interp
          18             :
          19             : contains
          20           0 :   subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight_distance)
          21             : !------------------------------------------------------------------------------------------------------------
          22             : ! This program computes weighting functions to map a variable of (im1,jm1) resolution to (im2,jm2) resolution
          23             : ! weight_x(im2,im1) is the weighting function for zonal interpolation
          24             : ! weight_y(jm2,jm1) is the weighting function for meridional interpolation
          25             : !
          26             : ! Author: Chih-Chieh (Jack) Chen  -- May 2010
          27             : !
          28             : !------------------------------------------------------------------------------------------------------------
          29             :   implicit none
          30             :   integer,  intent(in)  :: im1, jm1, im2, jm2
          31             :   logical,  intent(in)  :: use_flight_distance   !.true. = flight distance, .false. =  all mixing ratios
          32             :   real(r8), intent(in)  :: lon0(im1), lat0(jm1)
          33             :   real(r8), intent(out) :: weight_x(im2,im1), weight_y(jm2,jm1)
          34             :
          35           0 :   real(r8) :: lon1(im1), lat1(jm1)
          36           0 :   real(r8) :: lon2(im2), lat2(jm2)
          37           0 :   real(r8) :: slon1(im1+1), slon2(im2+1), slat1(jm1+1), slat2(jm2+1)
          38             :   real(r8) :: x1_west, x1_east, x2_west, x2_east
          39             :   real(r8) :: y1_south, y1_north, y2_south, y2_north
          40             :   integer  :: i1, j1, i2, j2
          41             :
          42           0 :   weight_x(:,:) = 0.0_r8
          43           0 :   weight_y(:,:) = 0.0_r8
          44             :
          45             : ! lon0 & lat0 are longitude & latitude on the source mesh in radians
          46             : ! convert lon1, lat1 from radians to degrees
          47           0 :   lon1(:) = lon0(:)/SHR_CONST_PI*180.0_r8
          48           0 :   lat1(:) = lat0(:)/SHR_CONST_PI*180.0_r8
          49             :
          50             : ! set up lon2, lat2 (target mesh), in CAM convention
          51           0 :   do i2=1,im2
          52           0 :    lon2(i2) = (float(i2)-1.0_r8)*360.0_r8/float(im2)
          53             :   enddo
          54           0 :   do j2=1,jm2
          55           0 :    lat2(j2) = -90.0_r8+(float(j2)-1.0_r8)*180.0_r8/(float(jm2)-1.0_r8)
          56             :   enddo
          57             :
          58             :
          59             : ! set up staggered longitudes (cell edges in x)
          60           0 :   do i1=2,im1
          61           0 :         slon1(i1) = (lon1(i1-1)+lon1(i1))/2.0_r8
          62             :   enddo
          63           0 :   slon1(1) = lon1(1)-(lon1(2)-lon1(1))/2.0_r8
          64           0 :   slon1(im1+1) = lon1(im1)+(lon1(im1)-lon1(im1-1))/2.0_r8
          65             :
          66           0 :   do i2=2,im2
          67           0 :         slon2(i2) = (lon2(i2-1)+lon2(i2))/2.0_r8
          68             :   enddo
          69           0 :   slon2(1) = lon2(1)-(lon2(2)-lon2(1))/2.0_r8
          70           0 :   slon2(im2+1) = lon2(im2)+(lon2(im2)-lon2(im2-1))/2.0_r8
          71             :
          72             : ! set up staggered lattiudes (cell edges in y)
          73           0 :   slat1(1)=-90.0_r8
          74           0 :   do j1=2,jm1
          75           0 :         slat1(j1) = (lat1(j1-1)+lat1(j1))/2.0_r8
          76             :   enddo
          77           0 :   slat1(jm1+1)=90.0_r8
          78             :
          79           0 :   slat2(1)=-90.0_r8
          80           0 :   do j2=2,jm2
          81           0 :         slat2(j2)=(lat2(j2-1)+lat2(j2))/2.0_r8
          82             :   enddo
          83           0 :   slat2(jm2+1)=90.0_r8
          84             :
          85             : ! compute Guassian weight for two meshes (discrete form of cos(lat).)
          86           0 :   do j1=1,jm1
          87           0 :         gw1(j1) = sin(slat1(j1+1)/180.0_r8*SHR_CONST_PI)-sin(slat1(j1)/180.0_r8*SHR_CONST_PI)
          88             :   enddo
          89             :
          90           0 :   do j2=1,jm2
          91           0 :         gw2(j2) = sin(slat2(j2+1)/180.0_r8*SHR_CONST_PI)-sin(slat2(j2)/180.0_r8*SHR_CONST_PI)
          92             :   enddo
          93             :
          94             :
          95             : ! add 360 to slon1 and slon2
          96           0 :   slon1(:) = slon1(:)+360.0_r8
          97           0 :   slon2(:) = slon2(:)+360.0_r8
          98             :
          99           0 :  do i2=1,im2
         100             :
         101             : ! target grid east-west boundaries
         102           0 :   x2_west=slon2(i2)
         103           0 :   x2_east=slon2(i2+1)
         104             :
         105           0 :   do i1=1,im1
         106             :
         107             : ! source grid east-west boundaries
         108           0 :    x1_west=slon1(i1)
         109           0 :    x1_east=slon1(i1+1)
         110             :
         111             : ! check if there is any overlap between the source grid and the target grid
         112             : ! if no overlap, then weighting is zero
         113             : ! there are three scenarios overlaps can take place
         114           0 :    if( (x1_west>=x2_west).and.(x1_east<=x2_east) ) then
         115             : ! case 1:
         116             : !                x1_west             x1_east
         117             : !                  |-------------------|
         118             : !            |---------------------------------|
         119             : !          x2_west                           x2_east
         120           0 :     if(use_flight_distance) then
         121           0 :      weight_x(i2,i1) = 1.0_r8
         122             :     else
         123           0 :      weight_x(i2,i1) =  (x1_east-x1_west)/(x2_east-x2_west)
         124             :     endif
         125           0 :    elseif ( (x1_west>=x2_west).and.(x1_west<x2_east) ) then
         126             : ! case 2:
         127             : !                x1_west                          x1_east
         128             : !                  |--------------------------------|
         129             : !            |---------------------------------|
         130             : !          x2_west                           x2_east
         131           0 :    if(use_flight_distance) then
         132           0 :      weight_x(i2,i1) = (x2_east-x1_west)/(x1_east-x1_west)
         133             :    else
         134           0 :      weight_x(i2,i1) = (x2_east-x1_west)/(x2_east-x2_west)
         135             :    endif
         136           0 :    elseif ( (x1_east>x2_west).and.(x1_east<=x2_east) ) then
         137             : ! case 3:
         138             : !       x1_west                          x1_east
         139             : !         |--------------------------------|
         140             : !                |---------------------------------|
         141             : !              x2_west                           x2_east
         142           0 :    if(use_flight_distance) then
         143           0 :      weight_x(i2,i1) = (x1_east-x2_west)/(x1_east-x1_west)
         144             :    else
         145           0 :      weight_x(i2,i1) = (x1_east-x2_west)/(x2_east-x2_west)
         146             :    endif
         147           0 :    elseif ( (x1_east>x2_east).and.(x1_west<x2_west) ) then
         148             : ! case 4:
         149             : !       x1_west                          x1_east
         150             : !         |--------------------------------|
         151             : !                |--------------------|
         152             : !              x2_west              x2_east
         153           0 :      weight_x(i2,i1) = 1._r8
         154             :    endif
         155             :
         156             :    enddo
         157             :   enddo
         158             :
         159             :
         160             : ! consider end points
         161           0 :       if(slon1(im1+1)>slon2(im2+1)) then
         162             : ! case 1:
         163             : !           slon1(im1)                slon1(im1+1) <--- end point
         164             : !              |-------------------------|
         165             : !           |----------------|......................|
         166             : !        slon2(im2)         slon2(im2+1)        slon2(2)  (note: slon2(im2+1) = slon2(1))
         167           0 :          if(use_flight_distance) then
         168           0 :             weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon1(im1+1)-slon1(im1))
         169             :          else
         170           0 :             weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon2(2)-slon2(1))
         171             :          endif
         172             :       endif
         173             :
         174           0 :       if(slon1(im1+1)<slon2(im2+1)) then
         175             : ! case 1:
         176             : !           slon1(im1)                slon1(im1+1)                  slon1(2)    (note: slon1(im1+1) = slon1(1))
         177             : !              |-------------------------|.............................|
         178             : !                   |-------------------------------|
         179             : !               slon2(im2)                        slon2(im2+1) <--- end point
         180           0 :          if(use_flight_distance) then
         181           0 :             weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon1(2)-slon1(1))
         182             :          else
         183           0 :             weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon2(2)-slon2(1))
         184             :          endif
         185             :       endif
         186             :
         187             :
         188             :
         189           0 :       do j2=1,jm2
         190             : ! target grid north-south boundaries
         191           0 :          y2_south=slat2(j2)
         192           0 :          y2_north=slat2(j2+1)
         193             :
         194           0 :          do j1=1,jm1
         195             :
         196             : ! source grid north-south boundaries
         197           0 :             y1_south=slat1(j1)
         198           0 :             y1_north=slat1(j1+1)
         199             :
         200             : ! check if there is any overlap between the source grid and the target grid
         201             : ! if no overlap, then weighting is zero
         202             : ! there are three scenarios overlaps can take place
         203             : ! note: there is Guassian weight to consider in the meridional direction!
         204             :
         205           0 :             if( (y1_south>=y2_south).and.(y1_north<=y2_north) ) then
         206             : ! case 1:
         207             : !                y1_south             y1_north
         208             : !                  |-------------------|
         209             : !            |---------------------------------|
         210             : !          y2_south                           y2_north
         211           0 :                 if(use_flight_distance) then
         212           0 :                    weight_y(j2,j1) =  1.0_r8
         213             :                 else
         214           0 :                    weight_y(j2,j1) =  gw1(j1)/gw2(j2)
         215             :                 endif
         216           0 :             elseif ( (y1_south>=y2_south).and.(y1_south<y2_north) ) then
         217             : ! case 2:
         218             : !                y1_south                          y1_north
         219             : !                  |--------------------------------|
         220             : !            |---------------------------------|
         221             : !          y2_south                           y2_north
         222           0 :                 if(use_flight_distance) then
         223           0 :                    weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south)
         224             :                 else
         225           0 :                    weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2)
         226             :                 endif
         227           0 :             elseif ( (y1_north>y2_south).and.(y1_north<=y2_north) ) then
         228             : ! case 3:
         229             : !       y1_south                          y1_north
         230             : !         |--------------------------------|
         231             : !                |---------------------------------|
         232             : !              y2_south                           y2_north
         233           0 :                 if(use_flight_distance) then
         234           0 :                    weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)
         235             :                 else
         236           0 :                    weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2)
         237             :                 endif
         238           0 :             elseif ( (y1_north>y2_north).and.(y1_south<y2_south) ) then
         239             : ! case 4:
         240             : !       y1_south                          y1_north
         241             : !         |--------------------------------|
         242             : !                |---------------------|
         243             : !              y2_south             y2_north
         244           0 :                 if(use_flight_distance) then
         245           0 :                    weight_y(j2,j1) = 1._r8
         246             :                 else
         247           0 :                    weight_y(j2,j1) = gw1(j1)/gw2(j2)
         248             :                 endif
         249             :             endif
         250             :
         251             :           enddo
         252             :         enddo
         253             :
         254           0 :   end subroutine xy_interp_init
         255             :
         256           0 :   subroutine xy_interp(im1,jm1,km1,im2,jm2,pcols,ncols,weight_x,weight_y,var_src,var_trg,lons,lats,count_x,count_y,index_x,index_y)
         257             : !-------------------------------------------------------------------------------------------------------------
         258             : ! This program interpolates var_src(im1,jm1,km1) to var_trg(im2,jm2,km1) based on weighting functions weight_x & weight_y.
         259             : !-------------------------------------------------------------------------------------------------------------
         260             :   implicit none
         261             :   integer,  intent(in)  :: im1   ! source number of longitudes
         262             :   integer,  intent(in)  :: jm1   ! source number of latitudes
         263             :   integer,  intent(in)  :: km1   ! source/target number of levels
         264             :   integer,  intent(in)  :: im2   ! target number of longitudes
         265             :   integer,  intent(in)  :: jm2   ! target number of latitudes
         266             :   integer,  intent(in)  :: pcols
         267             :   integer,  intent(in)  :: ncols
         268             :   real(r8), intent(in)  :: weight_x(im2,im1), weight_y(jm2,jm1)
         269             :   real(r8), intent(in)  :: var_src(im1,jm1,km1)
         270             :   integer,  intent(in)  :: lons(pcols), lats(pcols)
         271             :   integer,  intent(in)  :: count_x(im2), count_y(jm2)
         272             :   integer,  intent(in)  :: index_x(im2,im1), index_y(jm2,jm1)
         273             :   real(r8), intent(out) :: var_trg(pcols,km1)
         274             :   integer  :: n, i1, j1, k1, i2, j2, ii, jj
         275             :   real(r8) :: sum_x
         276             :
         277           0 :   var_trg(:,:) = 0.0_r8
         278             :
         279             :
         280           0 :   do k1=1,km1
         281           0 :      do n=1,ncols
         282             : ! interpolate in x
         283           0 :         do jj=1,count_y(lats(n))
         284           0 :            sum_x = 0.0_r8
         285           0 :            do ii=1,count_x(lons(n))
         286           0 :               sum_x = sum_x + var_src(index_x(lons(n),ii),index_y(lats(n),jj),k1)* &
         287           0 :                              weight_x(lons(n),index_x(lons(n),ii))
         288             :            enddo
         289           0 :            var_trg(n,k1) = var_trg(n,k1)+sum_x*weight_y(lats(n),index_y(lats(n),jj))
         290             :         enddo
         291             :      enddo
         292             :   enddo
         293             :
         294           0 :   end subroutine xy_interp
         295             :
         296             :
         297             : end module horizontal_interpolate

Generated by: LCOV version 1.14