! !################################################################## !################################################################## !###### ###### !###### 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,29 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 == 3) THEN ! HDF 4 format CALL hdfopen(hisfile(1:lenfil),1,nchr) CALL hdfrdi(nchr,"nx",nx,ireturn) CALL hdfrdi(nchr,"ny",ny,ireturn) CALL hdfrdi(nchr,"nz",nz,ireturn) CALL hdfrdi(nchr,"nzsoil",nzsoil,ireturn) CALL hdfrdi(nchr,"nstyp", nstyps,ireturn) CALL hdfrdi(nchr,"month", month, ireturn) CALL hdfrdi(nchr,"day", day, ireturn) CALL hdfrdi(nchr,"year", year, ireturn) CALL hdfrdi(nchr,"hour", hour, ireturn) CALL hdfrdi(nchr,"minute",minute,ireturn) CALL hdfrdi(nchr,"second",second,ireturn) CALL hdfrdi(nchr,"mapproj",mapproj,ireturn) CALL hdfrdr(nchr,"trulat1",trulat1,ireturn) CALL hdfrdr(nchr,"trulat2",trulat2,ireturn) CALL hdfrdr(nchr,"trulon", trulon, ireturn) CALL hdfrdr(nchr,"sclfct", sclfct, ireturn) CALL hdfrdr(nchr,"ctrlat", ctrlat, ireturn) CALL hdfrdr(nchr,"ctrlon", ctrlon, ireturn) CALL hdfrdr(nchr,"time",time,ireturn) CALL hdfrdr(nchr,"dx",dx,ireturn) CALL hdfrdr(nchr,"dy",dy,ireturn) CALL hdfclose(nchr,ireturn) 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, & rdummy,rdummy,rdummy,rdummy,rdummy, & 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 grid_id,parent_id,i_parent_start,j_parent_start, & i_parent_end,j_parent_end,parent_grid_ratio, & nx,ny,dx,dy, & mapproj,trulat1,trulat2,trulon,ctrlat_moad, & 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) :: grid_id, parent_id INTEGER, INTENT(IN) :: i_parent_start,j_parent_start INTEGER, INTENT(IN) :: i_parent_end, j_parent_end INTEGER, INTENT(IN) :: parent_grid_ratio INTEGER, INTENT(IN) :: nx,ny REAL, INTENT(IN) :: dx,dy INTEGER, INTENT(IN) :: mapproj REAL, INTENT(IN) :: trulat1,trulat2,trulon REAL, INTENT(IN) :: ctrlat_moad 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', grid_id, istatus) CALL put_dom_ti_integer(nfid,io_form, 'PARENT_ID', parent_id, istatus) CALL put_dom_ti_integer(nfid,io_form, 'I_PARENT_START',i_parent_start, istatus) CALL put_dom_ti_integer(nfid,io_form, 'J_PARENT_START',j_parent_start, istatus) CALL put_dom_ti_integer(nfid,io_form, 'I_PARENT_END', i_parent_end,istatus) CALL put_dom_ti_integer(nfid,io_form, 'J_PARENT_END', i_parent_end,istatus) CALL put_dom_ti_integer(nfid,io_form, 'PARENT_GRID_RATIO',parent_grid_ratio,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_moad,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*256 :: path_soiltemp_1deg INTEGER :: istatus CHARACTER(LEN=256) :: 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