!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE get_arps_dims_atts ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE get_arps_dims_atts(nchr,hinfmt,hisfile,lenfil, & 1,6
nx,ny,nz,nzsoil,nstyps, &
year,month,day,hour,minute,second,time, &
mapproj, sclfct,trulat1,trulat2,trulon, &
ctrlat,ctrlon,dx,dy,ireturn)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Get ARPS dimensions and attributes from ARPS time-independent
! history file.
!
!-----------------------------------------------------------------------
!
! Author: Yunheng Wang (09/01/2005)
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(OUT) :: nchr
INTEGER, INTENT(IN) :: hinfmt
CHARACTER(*), INTENT(IN) :: hisfile
INTEGER, INTENT(IN) :: lenfil
INTEGER, INTENT(OUT) :: nx,ny,nz,nzsoil,nstyps
INTEGER, INTENT(OUT) :: year, month,day,hour,minute,second
REAL, INTENT(OUT) :: time
INTEGER, INTENT(OUT) :: mapproj
REAL, INTENT(OUT) :: sclfct
REAL, INTENT(OUT) :: trulat1,trulat2,trulon
REAL, INTENT(OUT) :: ctrlat,ctrlon
REAL, INTENT(OUT) :: dx,dy
INTEGER, INTENT(OUT) :: ireturn
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
INTEGER :: idummy
REAL :: rdummy
CHARACTER(LEN=40) :: fmtverin
CHARACTER(LEN=10) :: tmunit
CHARACTER(LEN=80) :: runname
CHARACTER(LEN=12) :: label
INTEGER :: nocmnt, totin
CHARACTER(LEN=80) :: cmnt(50)
REAL :: thisdmp, tstop
REAL :: latitud, xgrdorg, ygrdorg, umove, vmove
REAL, ALLOCATABLE :: x(:)
REAL, ALLOCATABLE :: y(:)
INTEGER :: i
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Begin of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF( hinfmt == 1 ) THEN
CALL getunit
( nchr )
OPEN(UNIT=nchr,FILE=hisfile(1:lenfil),STATUS='old', &
FORM='unformatted',IOSTAT=ireturn)
READ(nchr) fmtverin
READ(nchr) runname
READ(nchr) nocmnt
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
READ(nchr) cmnt(i)
END DO
END IF
READ(nchr) time,tmunit
READ(nchr) nx, ny, nz,nzsoil
READ(nchr) idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,totin, &
idummy,idummy,idummy,mapproj,month, &
day, year, hour, minute, second
READ(nchr) umove, vmove, xgrdorg,ygrdorg,trulat1, &
trulat2,trulon,sclfct,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
tstop,thisdmp,latitud,ctrlat,ctrlon
IF ( totin /= 0 ) THEN
READ(nchr) nstyps,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
IF ( nstyps < 1 ) nstyps = 1
READ(nchr) rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
END IF
ALLOCATE(x(nx), STAT = idummy)
ALLOCATE(y(ny), STAT = idummy)
READ(nchr) label
READ(nchr) x
READ(nchr) label
READ(nchr) y
CLOSE (nchr)
CALL retunit( nchr )
dx = x(2) - x(1)
dy = y(2) - y(1)
DEALLOCATE(x,y)
ELSE IF (hinfmt == 7 .OR. hinfmt == 8) THEN ! NetCDF format
CALL netopen(hisfile(1:lenfil),'R',nchr)
CALL net_getdims(nchr,nx,ny,nz,nzsoil,nstyps,ireturn)
CALL net_getatts(nchr,runname,nocmnt,cmnt,dx,dy, &
year,month,day,hour,minute,second,thisdmp,tstop, &
mapproj,sclfct,trulat1,trulat2,trulon,latitud, &
ctrlat,ctrlon,xgrdorg,ygrdorg,umove,vmove, &
idummy,idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,ireturn)
! CALL netreadTime(nchr,1,'Time',time)
time = 0.0
CALL netclose(nchr)
ELSE
WRITE(6,'(a,i3,a)') &
' Data format flag had an invalid value ', &
hinfmt ,' program stopped.'
CALL arpsstop('arpsstop called from get_arps_dims_atts wrong flag',1)
END IF
RETURN
END SUBROUTINE get_arps_dims_atts
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE write_static_attribute ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE write_static_attribute(nfid,io_form,start_date, & 1,32
nx,ny,dx,dy, &
mapproj,trulat1,trulat2,trulon,ctrlat,ctrlon, &
lat_ll,lat_ul,lat_ur,lat_lr, &
lon_ll,lon_ul,lon_ur,lon_lr, istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Write attributes of WRF static file.
!
!-----------------------------------------------------------------------
USE wrf_metadata
IMPLICIT NONE
INTEGER, INTENT(IN) :: nfid
INTEGER, INTENT(IN) :: io_form
CHARACTER(*), INTENT(IN) :: start_date
INTEGER, INTENT(IN) :: nx,ny
REAL, INTENT(IN) :: dx,dy
INTEGER, INTENT(IN) :: mapproj
REAL, INTENT(IN) :: trulat1,trulat2,trulon
REAL, INTENT(IN) :: ctrlat, ctrlon
REAL, INTENT(IN) :: lat_ll(4),lat_ul(4),lat_ur(4),lat_lr(4)
REAL, INTENT(IN) :: lon_ll(4),lon_ul(4),lon_ur(4),lon_lr(4)
INTEGER, INTENT(OUT) :: istatus
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
INTEGER :: map_proj
CHARACTER(LEN=20) :: map_string
REAL :: latcorns(16), loncorns(16)
INTEGER :: i,k
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Begin of executable code ...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
map_string(1:20) = ' '
IF(ABS(mapproj) == 1) THEN
map_proj = 2
map_string = 'POLAR STEREOGRAPHIC'
ELSE IF(ABS(mapproj) == 2) THEN
map_proj = 1
map_string = 'LAMBERT CONFORMAL'
ELSE IF(ABS(mapproj) == 3) THEN
map_proj = 3
map_string = 'MERCATOR'
ELSE
WRITE(6,*) 'Unknown map projection, ', mapproj
STOP
END IF
DO i = 1,4
k = (i-1)*4+1
latcorns(k) = lat_ll(i)
loncorns(k) = lon_ll(i)
k = k+1
latcorns(k) = lat_ul(i)
loncorns(k) = lon_ul(i)
k = k+1
latcorns(k) = lat_ur(i)
loncorns(k) = lon_ur(i)
k = k+1
latcorns(k) = lat_lr(i)
loncorns(k) = lon_lr(i)
END DO
IF (io_form == 7) CALL enter_ncd_define(nfid,istatus)
CALL put_dom_ti_char(nfid,io_form,'TITLE', &
'WRF static file from ARPS2WRF',istatus)
CALL put_dom_ti_char(nfid,io_form,'simulation_name', &
'WRFSTATIC',istatus)
CALL put_dom_ti_char(nfid,io_form,'START_DATE', &
start_date(1:19),istatus)
CALL put_dom_ti_integer(nfid,io_form,'DYN_OPT', 2,istatus)
CALL put_dom_ti_integer(nfid,io_form,'WEST-EAST_GRID_DIMENSION', &
nx, istatus)
CALL put_dom_ti_integer(nfid,io_form,'SOUTH-NORTH_GRID_DIMENSION', &
ny, istatus)
CALL put_dom_ti_integer(nfid,io_form, 'GRID_ID', 1, istatus)
CALL put_dom_ti_integer(nfid,io_form, 'PARENT_ID', 1, istatus)
CALL put_dom_ti_integer(nfid,io_form, 'I_PARENT_START',1, istatus)
CALL put_dom_ti_integer(nfid,io_form, 'J_PARENT_START',1, istatus)
CALL put_dom_ti_integer(nfid,io_form, 'I_PARENT_END', nx,istatus)
CALL put_dom_ti_integer(nfid,io_form, 'J_PARENT_END', ny,istatus)
CALL put_dom_ti_integer(nfid,io_form, 'PARENT_GRID_RATIO',1,istatus)
CALL put_dom_ti_real (nfid,io_form, 'DX',dx,istatus)
CALL put_dom_ti_real (nfid,io_form, 'DY',dy,istatus)
CALL put_dom_ti_varreal(nfid,io_form, 'corner_lats',latcorns,16,istatus)
CALL put_dom_ti_varreal(nfid,io_form, 'corner_lons',loncorns,16,istatus)
CALL put_dom_ti_real (nfid,io_form, 'CEN_LAT',ctrlat,istatus)
CALL put_dom_ti_real (nfid,io_form, 'CEN_LON',ctrlon,istatus)
CALL put_dom_ti_integer(nfid,io_form, 'FLAG_STATIC',1,istatus)
CALL put_dom_ti_char(nfid,io_form,'map_projection',TRIM(map_string), &
istatus)
CALL put_dom_ti_integer(nfid,io_form, 'MAP_PROJ', map_proj,istatus)
CALL put_dom_ti_real (nfid,io_form, 'MOAD_CEN_LAT', ctrlat,istatus)
CALL put_dom_ti_real (nfid,io_form, 'STAND_LON', trulon,istatus)
CALL put_dom_ti_real (nfid,io_form, 'TRUELAT1', trulat1,istatus)
CALL put_dom_ti_real (nfid,io_form, 'TRUELAT2', trulat2,istatus)
CALL put_dom_ti_integer(nfid,io_form, 'ISWATER',ISWATER,istatus)
CALL put_dom_ti_integer(nfid,io_form, 'ISICE', ISICE, istatus)
CALL put_dom_ti_char (nfid,io_form, 'MMINLU', 'USGS', istatus)
IF (io_form == 7) CALL exit_ncd_define(nfid,istatus)
RETURN
END SUBROUTINE write_static_attribute
!
! ================================================================
!
SUBROUTINE RDGEODAT(n2,n3,xt,yt,deltax,deltay, & 5,1
ofn,wvln,silwt,maxdatacat, &
datr,dats,datln,datlt,which_data, istat)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n2,n3
REAL, INTENT(IN) :: deltax,deltay,wvln,silwt
! wvln = TOPTWVL_PARM_WRF from wrfsi.nl section &hgridspec
! siwt = SILAVWT_PARM_WRF
REAL, INTENT(IN) :: xt(n2),yt(n3)
INTEGER, INTENT(IN) :: maxdatacat
CHARACTER(LEN=*),INTENT(IN) :: OFN
REAL, INTENT(OUT):: datr(n2,n3)
REAL, INTENT(OUT):: dats(n2,n3,maxdatacat)
REAL, INTENT(OUT):: datln(n2,n3)
REAL, INTENT(OUT):: datlt(n2,n3)
LOGICAL, INTENT(OUT):: which_data
INTEGER, INTENT(OUT):: istat
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
INTEGER, PARAMETER :: IODIM = 5800000
! this parameter can be increased if we ever read data with
! resolution finer than 30 sec, or if the tilesize for 30s
! data becomes greater than 10x10 deg.
!
INTEGER :: lb,lbf
INTEGER :: mof,np,niq,njq,nx,ny
REAL :: deltallo,deltaxq,deltayq, deltaxp,deltayp
INTEGER :: no, iblksizo, isbego, iwbego
REAL :: rsoff,rwoff
INTEGER :: lcat
CHARACTER(LEN=180) ::TITLE
REAL :: std_lon = -100.0 ! not used anywhere
REAL, PARAMETER :: erad = 6371200.0 ! not used anywhere
CHARACTER(LEN=180) :: tmpstr
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! begin of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! *********************
nx = n2-1
ny = n3-1
! **********************
lcat = 1
lbf = LEN_TRIM(ofn)
tmpstr(1:180) = ' '
write(tmpstr,'(a)') ofn(1:lbf)
title = ofn(1:lbf)//'HEADER'
lb = LEN_TRIM(title)
OPEN(29,STATUS='OLD',FILE=title(1:lb),FORM='FORMATTED')
READ(29,*) iblksizo,no,isbego,iwbego,rsoff,rwoff
print *,'title = ',title
print *,'rsoff,rwoff = ',rsoff,rwoff
print *,'isbego,iwbego = ',isbego,iwbego
print *,'iblksizo,no = ',iblksizo,no
CLOSE(29)
IF (no > 1) THEN
deltallo = FLOAT(iblksizo)/FLOAT(no-1)
ELSE IF(no == 1) THEN
deltallo = FLOAT(iblksizo)/FLOAT(no)
ELSE
print*,'HEADER value NO = 0'
istat = -1
RETURN
ENDIF
mof = iodim/(no*no)
IF (ofn(lbf:lbf) == 'G' .OR. ofn(lbf:lbf) == 'A') THEN ! Green_frac or Albedo
lcat=12
IF (no == 1250) mof = 1
END IF
! MOF determines the number of files held in buffer while reading
! DEM data; it saves some time when buffer data can be used instead
! of reading DEM file again. Originally MOF was 4.
IF (mof > 10) mof = 5
deltaxq = 0.5*wvln*deltax
deltayq = 0.5*wvln*deltay
print *,'deltaxq,deltayq=',deltaxq,deltayq
niq = INT(FLOAT(nx)*deltax/deltaxq) + 4
njq = INT(FLOAT(ny)*deltay/deltayq) + 4
np = MIN(10,MAX(1,INT(deltaxq/(deltallo*111000.))))
print *,' np=',np
deltaxp = deltaxq/FLOAT(np)
deltayp = deltayq/FLOAT(np)
CALL SFCOPQR(NO,MOF,NP,NIQ,NJQ,N2,N3,lcat, &
XT,YT,90.,std_lon,ERAD,RWOFF,RSOFF, &
DELTALLO,DELTAXP,DELTAYP,DELTAXQ,DELTAYQ, &
IBLKSIZO,ISBEGO,IWBEGO,DATR,DATS,DATLN,DATLT, &
tmpstr,WVLN,SILWT,which_data,maxdatacat,istat)
RETURN
END SUBROUTINE rdgeodat
SUBROUTINE get_path_to_soiltemp_1deg (path_soiltemp_1deg,istatus) 1
IMPLICIT NONE
CHARACTER*200 :: path_soiltemp_1deg
INTEGER :: istatus
CHARACTER(LEN=200) :: path_to_soiltemp_1deg
COMMON /lapscmn/ path_to_soiltemp_1deg
path_soiltemp_1deg = path_to_soiltemp_1deg
RETURN
END SUBROUTINE
SUBROUTINE interpolate_masked_val(nx_src, ny_src, lmask_src, data_src, & 1
nx_out, ny_out, lmask_out, data_out, &
isrcr, jsrcr, make_srcmask, &
min_val, max_val, def_val, val_mask, &
method)
IMPLICIT none
INTEGER, INTENT(IN) :: nx_src
INTEGER, INTENT(IN) :: ny_src
REAL, INTENT(INOUT) :: lmask_src(nx_src,ny_src)
REAL, INTENT(IN) :: data_src(nx_src,ny_src)
INTEGER, INTENT(IN) :: nx_out, ny_out
REAL, INTENT(IN) :: lmask_out(nx_out,ny_out)
REAL, INTENT(INOUT) :: data_out(nx_out,ny_out)
REAL, INTENT(IN) :: isrcr(nx_out,ny_out)
REAL, INTENT(IN) :: jsrcr(nx_out,ny_out)
LOGICAL, INTENT(IN) :: make_srcmask
REAL, INTENT(IN) :: min_val, max_val, def_val, val_mask
INTEGER, INTENT(IN) :: method
INTEGER :: i,j,k,l,ilo,jlo, ii, jj, kk, ll
INTEGER :: search_rad
INTEGER :: search_rad_max
REAL, PARAMETER :: bad_point = -99999.
LOGICAL :: bad_points_flag
LOGICAL :: fixed_bad
INTEGER :: total_points
INTEGER :: points_val_mask
INTEGER :: points_skipped
INTEGER :: num_to_fix
INTEGER :: fixed_with_out
INTEGER :: fixed_with_src
INTEGER :: fixed_with_def
INTEGER :: close_pts
INTEGER :: ic,jc,iic,jjc
REAL :: deltagx, deltagy
REAL :: dix, djy
REAL :: data_close(4)
REAL :: distance(4)
REAL :: distx,disty
REAL :: sum_dist
REAL :: a,b,c,d,e,f,g,h
REAL :: stl(4,4)
REAL :: valb
REAL :: srcmask(nx_src,ny_src)
LOGICAL :: wrapped_source
INTEGER, PARAMETER :: METHOD_NEAREST = 0
INTEGER, PARAMETER :: METHOD_LINEAR = 1
INTEGER, PARAMETER :: METHOD_HIGHER = 2
REAL :: oned
IF ((MINVAL(isrcr).LT.1.0).OR. (MAXVAL(isrcr).GT.nx_src))THEN
wrapped_source = .TRUE.
ELSE
wrapped_source = .FALSE.
ENDIF
! Can we use the source data land mask that was input, or do we need to make it
! from the min_val/max_val arguments?
IF (make_srcmask) THEN
PRINT '(A)', 'INTERPOLATE_MASKED_VAL: Making source data land mask...'
! This code has the flexibility to handled either water data points (lmask =0)
! or land points (lmask = 1). So if we are making the landmask using valid
! range, but we do not know anything about the variable other than the
! valid mask value, we need to figure out which points are 0 and which should be
! 1. To do this, initialize the array to the invalid value, which we determine
! from the valid value.
IF (val_mask .EQ. 0) THEN
lmask_src(:,:) = 1
ELSEIF (val_mask .EQ. 1) THEN
lmask_src(:,:) = 0
ELSE
lmask_src(:,:) = -1
ENDIF
! Now figure out where to set the mask using a WHERE statement
! NOTE: In this next line, if val_mask is 2, then the lmask_src
! is going to be set to 2, so we need to be careful in some of the
! subsequent IF statements where lmask_src is compared to lmask_out
WHERE((data_src .GE. min_val).AND.(data_src .LE. max_val)) lmask_src = val_mask
ELSE
PRINT '(A)', 'INTERPOLATE_MASKED: Using source landmask field.'
ENDIF
bad_points_flag = .false.
! Initialize counters
total_points = 0
points_val_mask = 0
points_skipped = 0
num_to_fix = 0
fixed_with_out = 0
fixed_with_src = 0
fixed_with_def = 0
! Select interpolation method. Putting the case statement out here
! increases the amount of replicated code but should be more efficient
! than checking this condition at every grid point.
SELECT CASE(method)
CASE(METHOD_NEAREST)
! Use nearest neigbor
PRINT '(A)', 'Masked interpolation using nearest neighbor value...'
out_j_loop_1: DO j = 1, ny_out
out_i_loop_1: DO i = 1, nx_out
total_points = total_points + 1
! We only need to process this point if the lmask_out is equal
! to val_mask. For example, if one is processing soil parameters,
! which are only valid at land points (lmask = 1), then the user
! passes in 1. for the val_mask. Any output point that is water
! (lmask = 0) will then be skipped. During this loop, if we have
! points that cannot be properly assigned, we will mark them as bad
! and set the bad_points_lag. Exception is if val_mask = 2, which
! implies that we are doing masked interpolation for both land
! and water, but allowing only water points to influence water
! points and land points to influence landpoints.
IF ((lmask_out(i,j) .EQ. val_mask).OR.(val_mask .EQ. 2)) THEN
! Process this point
points_val_mask = points_val_mask + 1
ilo = NINT(isrcr(i,j))
! Account for horizontal wrapping as in the case
! of global data sets. This is the only time
! it is possible for ilo < 1 or ilo > nx_src
IF (ilo .LT.1) ilo = nx_src
IF (ilo .GT. nx_src) ilo = 1
jlo = NINT(jsrcr(i,j))
! See if this point can be used
IF ((lmask_src(ilo,jlo).EQ. lmask_out(i,j)).OR. &
(lmask_src(ilo,jlo).EQ.2)) THEN
data_out(i,j) = data_src(ilo,jlo)
ELSE
data_out(i,j) = bad_point
num_to_fix = num_to_fix + 1
bad_points_flag = .true.
ENDIF
ELSE
! The output grid does not require a value for this point
! But do not zero out in case this is a field begin
! done twice (once for water and once for land, e.g.
! SKINTEMP
!data_out(i,j) = 0.
points_skipped = points_skipped + 1
ENDIF
ENDDO out_i_loop_1
ENDDO out_j_loop_1
CASE (METHOD_LINEAR)
! Use a 4-point interpolation
PRINT '(A)', 'Masked interpolation using 4-pt linear interpolation...'
deltagx = 1.
deltagy = 1.
out_j_loop_2: DO j = 1, ny_out
out_i_loop_2: DO i = 1, nx_out
total_points = total_points + 1
! We only need to process this point if the lmask_out is equal
! to val_mask. For example, if one is processing soil parameters,
! which are only valid at land points (lmask = 1), then the user
! passes in 1. for the val_mask. Any output point that is water
! (lmask = 0) will then be skipped. During this loop, if we have
! points that cannot be properly assigned, we will mark them as bad
! and set the bad_points_lag.
IF ((lmask_out(i,j) .EQ. val_mask).OR.(val_mask .EQ. 2)) THEN
! Process this point
points_val_mask = points_val_mask + 1
! If ilo < 1 or > nx_src, this is a wrapped source dataset
ilo = FLOOR(isrcr(i,j))
IF (ilo .EQ. 0) ilo = nx_src
jlo = MIN(FLOOR(jsrcr(i,j)),ny_src-1)
dix = isrcr(i,j) - FLOAT(ilo)
IF (dix .LT.0) dix = dix + FLOAT(nx_src)
djy = jsrcr(i,j) - FLOAT(jlo)
! Loop around the four surrounding points
! and count up the number of points we can use based
! on common mask value
close_pts = 0
sum_dist = 0.
outer_four_j: DO jc = 0,1
outer_four_i: DO ic = 0,1
iic = ilo + ic
IF(iic .GT. nx_src) iic = iic - nx_src
jjc = jlo + jc
IF ((lmask_src(iic,jjc).EQ. lmask_out(i,j)).OR. &
(lmask_src(iic,jjc).EQ.2)) THEN
close_pts = close_pts + 1
data_close(close_pts) = data_src(iic,jjc)
! Compute distance to this valid point
! in grid units and add to sum of distances (actually,
! we are doing a weight, which is inversely proportional
! to distance)
IF (ic .EQ. 0) THEN
distx = deltagx - dix
ELSE
distx = dix
ENDIF
IF (jc .EQ. 0) THEN
disty = deltagy - djy
ELSE
disty = djy
ENDIF
distance(close_pts) = SQRT(distx**2+disty**2)
sum_dist = sum_dist + distance(close_pts)
ENDIF
ENDDO outer_four_i
ENDDO outer_four_j
! Did we find at least one point in the surrounding four
! that was usable?
IF (close_pts .GT. 0) THEN
! If we have all four points, then do bilinear interpolation
IF (close_pts .EQ. 4) THEN
data_out(i,j) = ((deltagx - dix)*((deltagy-djy)*data_close(1) &
+ djy*data_close(3)) &
+ dix*((deltagy-djy)*data_close(2) &
+ djy*data_close(4))) &
/ (deltagx * deltagy)
ELSE IF ((close_pts .GT. 1).AND.(close_pts .LT. 4)) THEN
! Simple distance-weighted average by computing
! the sum of all distances to each point and using
! each individual distance divided by the total
! distance as the weighting
data_out(i,j) = 0.
DO k = 1, close_pts
data_out(i,j) = data_out(i,j) + &
(distance(k)/sum_dist) * data_close(k)
ENDDO
ELSE
! Set output value = to one point we found
data_out(i,j) = data_close(1)
ENDIF
ELSE
bad_points_flag = .true.
data_out(i,j) = bad_point
num_to_fix = num_to_fix + 1
ENDIF
ELSE
! The output grid does not require a value for this point
! But do not zero out in case this is a field begin
! done twice (once for water and once for land, e.g.
! SKINTEMP
points_skipped = points_skipped + 1
ENDIF
ENDDO out_i_loop_2
ENDDO out_j_loop_2
CASE (METHOD_HIGHER)
! 16-point interpolation
out_j_loop_3: DO j = 1, ny_out
out_i_loop_3: DO i = 1, nx_out
total_points = total_points + 1
! We only need to process this point if the lmask_out is equal
! to val_mask. For example, if one is processing soil parameters,
! which are only valid at land points (lmask = 1), then the user
! passes in 1. for the val_mask. Any output point that is water
! (lmask = 0) will then be skipped. During this loop, if we have
! points that cannot be properly assigned, we will mark them as bad
! and set the bad_points_lag.
IF ((lmask_out(i,j) .EQ. val_mask).OR.(val_mask .EQ. 2)) THEN
! Process this point
points_val_mask = points_val_mask + 1
! Do a 4x4 loop around in the input data around the output
! point to get our 16 points of influence. Only use those
! that are appropriately masked
valb = 0.
close_pts = 0
ilo = INT(isrcr(i,j)+0.00001)
jlo = INT(jsrcr(i,j)+0.00001)
dix = isrcr(i,j) - ilo
djy = jsrcr(i,j) - jlo
IF ( (ABS(dix).GT.0.0001).OR.(ABS(djy).GT.0.0001) ) THEN
! Do the interpolation loop
stl(:,:) = 0.
loop_16_1: DO k = 1,4
kk = ilo + k - 2
IF ((kk .LT. 1).OR.(kk .GT. nx_src)) THEN
IF (.NOT. wrapped_source) THEN
CYCLE loop_16_1
ELSE
IF (kk .LT. 1) THEN
kk = kk + nx_src
ELSE
kk = kk - nx_src
ENDIF
ENDIF
ENDIF
loop_16_2: DO l = 1, 4
ll = jlo + l - 2
IF ((ll .LT. 1).OR.(ll .GT. ny_src)) CYCLE loop_16_2
! Check land mask at this source point
IF ((lmask_src(kk,ll).NE. lmask_out(i,j)).AND. &
(lmask_src(kk,ll).NE.2)) CYCLE loop_16_2
! If we are here, then mask tests passed
stl(k,l) = data_src(kk,ll)
IF ( (stl(k,l) .EQ. 0.).AND.(min_val.LE.0.).AND. &
(max_val.GE.0.) ) THEN
stl = 1.E-5
ENDIF
close_pts = close_pts + 1
ENDDO loop_16_2
ENDDO loop_16_1
! Did we find any valid points?
IF ( (close_pts .GT. 0).AND. ( &
(stl(2,2).GT.0.).AND.(stl(2,3).GT.0.).AND. &
(stl(3,2).GT.0.).AND.(stl(3,3).GT.0.) ) ) THEN
a = oned(dix,stl(1,1),stl(2,1),stl(3,1),stl(4,1))
b = oned(dix,stl(1,2),stl(2,2),stl(3,2),stl(4,2))
c = oned(dix,stl(1,3),stl(2,3),stl(3,3),stl(4,3))
d = oned(dix,stl(1,4),stl(2,4),stl(3,4),stl(4,4))
valb = oned(djy,a,b,c,d)
IF (close_pts .NE. 16) THEN
e = oned(djy,stl(1,1),stl(1,2),stl(1,3),stl(1,4))
f = oned(djy,stl(2,1),stl(2,2),stl(2,3),stl(2,4))
g = oned(djy,stl(3,1),stl(3,2),stl(3,3),stl(3,4))
h = oned(djy,stl(4,1),stl(4,2),stl(4,3),stl(4,4))
valb = (valb+oned(dix,e,f,g,h)) * 0.5
ENDIF
data_out(i,j) = valb
ELSE
bad_points_flag = .true.
data_out(i,j) = bad_point
num_to_fix = num_to_fix + 1
ENDIF
ELSE
! We are right on a source point, so try to use it
IF ((lmask_src(ilo,jlo).EQ.lmask_out(i,j)).OR.&
(lmask_src(ilo,jlo).EQ.2)) THEN
data_out(i,j) = data_src(ilo,jlo)
ELSE
bad_points_flag = .true.
data_out(i,j) = bad_point
num_to_fix = num_to_fix + 1
ENDIF
ENDIF
ELSE
! The output grid does not require a value for this point
!data_out(i,j) = 0.
points_skipped = points_skipped + 1
ENDIF
ENDDO out_i_loop_3
ENDDO out_j_loop_3
END SELECT
! Do we need to correct bad points?
IF (bad_points_flag) THEN
search_rad_max = 10
fix_bad_j: DO j = 1, ny_out
fix_bad_i: DO i = 1, nx_out
IF (data_out(i,j).NE. bad_point) CYCLE fix_bad_i
! First, search for nearest non-bad point in the output domain
! which is usually higher resolution.
fixed_bad = .false.
search_out_loop: DO search_rad = 1, search_rad_max
search_out_j: DO ll = -(search_rad-1), (search_rad-1),1
jj = j + ll
IF ((jj .LT. 1).OR.(jj .GT. ny_out)) CYCLE search_out_j
search_out_i: DO kk = -(search_rad), search_rad, 1
ii = i + kk
IF ((ii .LT. 1).OR.(ii .GT. nx_out)) CYCLE search_out_i
IF ((data_out(ii,jj).NE.bad_point).AND. &
(lmask_out(ii,jj) .EQ. lmask_out(i,j)) ) THEN
data_out(i,j) = data_out(ii,jj)
fixed_bad = .true.
fixed_with_out = fixed_with_out + 1
EXIT search_out_loop
ENDIF
ENDDO search_out_i
ENDDO search_out_j
ENDDO search_out_loop
! Did we fix the point? If not, then do same search on src data.
IF (.NOT. fixed_bad) THEN
search_rad_max = 10
search_src_loop: DO search_rad = 1, search_rad_max
search_src_j: DO ll = -(search_rad-1), (search_rad-1),1
jj = NINT(jsrcr(i,j)) + ll
IF ((jj .LT. 1).OR.(jj .GT. ny_src)) CYCLE search_src_j
search_src_i: DO kk = -(search_rad), search_rad, 1
ii = NINT(isrcr(i,j)) + kk
IF ((ii .LT. 1).OR.(ii .GT. nx_src)) THEN
IF (.NOT. wrapped_source) THEN
CYCLE search_src_i
ELSE
IF (ii.LT.1) THEN
ii = nx_src + ii
ELSE
ii = ii - nx_src
ENDIF
ENDIF
ENDIF
IF ((lmask_src(ii,jj).EQ.lmask_out(i,j)).OR. &
(lmask_src(ii,jj).EQ.2)) THEN
data_out(i,j) = data_src(ii,jj)
fixed_bad = .true.
fixed_with_src = fixed_with_src + 1
EXIT search_src_loop
ENDIF
ENDDO search_src_i
ENDDO search_src_j
ENDDO search_src_loop
ENDIF
! Now is the point fixed? If not, we have to use a default value.
IF (.NOT.fixed_bad) THEN
fixed_with_def = fixed_with_def + 1
data_out(i,j) = def_val
PRINT '(A,F10.3,A,2I5)', 'INTERPOLATE_MASKED: Bogus value of ', def_val, &
' used at point ', i, j
ENDIF
ENDDO fix_bad_i
ENDDO fix_bad_j
ENDIF
PRINT '(A)', '----------------------------------------'
PRINT '(A)', 'MASKED INTERPOLATION SUMMARY: '
PRINT '(A,I10)', ' TOTAL POINTS IN GRID: ', total_points
PRINT '(A,I10)', ' POINTS NEEDING VALUES: ', points_val_mask
PRINT '(A,I10)', ' POINTS NOT REQUIRED: ', points_skipped
PRINT '(A,I10)', ' POINTS NEEDING FIX: ', num_to_fix
PRINT '(A,I10)', ' POINTS FIXED WITH OUT GRID: ', fixed_with_out
PRINT '(A,I10)', ' POINTS FIXED WITH SRC GRID: ', fixed_with_src
PRINT '(A,I10)', ' POINTS FIXED WITH DEF VAL: ', fixed_with_def
PRINT '(A)', '----------------------------------------'
RETURN
END SUBROUTINE interpolate_masked_val
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
FUNCTION oned(x,a,b,c,d)
IMPLICIT NONE
REAL :: x,a,b,c,d,oned
oned = 0.
IF ( x .EQ. 0. ) THEN
oned = b
ELSE IF ( x .EQ. 1. ) THEN
oned = c
END IF
IF(b*c.NE.0.) THEN
IF ( a*d .EQ. 0. ) THEN
IF ( ( a .EQ. 0 ) .AND. ( d .EQ. 0 ) ) THEN
oned = b*(1.0-x)+c*x
ELSE IF ( a .NE. 0. ) THEN
oned = b+x*(0.5*(c-a)+x*(0.5*(c+a)-b))
ELSE IF ( d .NE. 0. ) THEN
oned = c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c))
END IF
ELSE
oned = (1.0-x)*(b+x*(0.5*(c-a)+x*(0.5*(c+a)-b)))+ &
x*(c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c)))
END IF
END IF
END FUNCTION oned
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE latlon_2_llij(np,glat,glon,lli,llj,lat0,lon0,dlat,dlon,cgrddef) 1
IMPLICIT NONE
INTEGER, INTENT(IN) :: np
REAL, INTENT(IN) :: glat(np), glon(np)
REAL, INTENT(OUT) :: lli(np), llj(np)
REAL, INTENT(IN) :: lat0, lon0
REAL, INTENT(IN) :: dlat, dlon
CHARACTER(LEN=1), INTENT(IN) :: cgrddef
INTEGER :: n
REAL :: diff
REAL :: dlond, dlatd
dlond=dlon
dlatd=dlat
IF (cgrddef .EQ. 'S') THEN
DO n=1,np
diff=glon(n)-lon0
IF (diff < 0.) diff=diff+360.
IF (diff >= 360.) diff=diff-360.
lli(n) = diff/dlond+1.
llj(n) = (glat(n)-lat0)/dlatd+1.
END DO
ELSE IF (cgrddef .EQ. 'N') THEN
DO n=1,np
diff=glon(n)-lon0
IF (diff < 0.) diff=diff+360.
IF (diff >= 360.) diff=diff-360.
lli(n) = diff/dlon+1.
llj(n) = (lat0-glat(n))/dlat+1.
END DO
ELSE
print*, 'you must specify whether the standard &
& lat is Southern or Northern boundary'
END IF
RETURN
END SUBROUTINE latlon_2_llij
FUNCTION bilinear_interp(i,j,imax,jmax,array_2d) RESULT (zo)
IMPLICIT NONE
INTEGER, INTENT(IN) :: i,j,imax,jmax
REAL, INTENT(IN) :: array_2d(imax,jmax)
REAL :: zo
REAL :: z1,z2,z3,z4
REAL :: fraci,fracj
fraci = 0.5
fracj = 0.5
z1 = array_2d(i , j )
z2 = array_2d(i-1, j )
z3 = array_2d(i-1, j-1)
z4 = array_2d(i , j-1)
zo = Z1 + (Z2-Z1)*fraci + (Z4-Z1)*fracj - (Z2+Z4-Z3-Z1)*fraci*fracj
END FUNCTION bilinear_interp