! !################################################################## !################################################################## !###### ###### !###### 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