----------------------------------------------------------------------- --- CALPUFF Modeling System --- ----------------------------------------------------------------------- Model Change Bulletin : G Dated : April 21, 2011 ----------------------------------------------------------------------- CALMET (From Version 5.812 Level 080512 to Version 5.815, Level 110421) ----------------------------------------------------------------------- -- Problem Area 1 -- When using multiple multi-hourly M3D files with overlapping times CALMET simulation stopped sooner than the requested end time. The fix does not affect results. Modified: RDMM5 -- Problem Area 2 -- When using convective overwater mixing height options (IMIXH>0) the overwater convective mixing height grows from the previous hour mixing height, not from the previous hour convective mixing height. This affects the values of overwater convective mixing heights overwater if conditions switch from stable/neutral to convective overwater during the simulation. Modified: WATER,WATERP -- Problem Area 3 -- At night when CALMET temperatures are computed from 3D.DAT temperatures (ITPROG>0) and CALMET lower level(s) are below the lowest 3D.DAT level, CALMET lowest temperature (s), which in that case are based on vertical extrapolation of the lowest 3D.DAT temperatures, assumed radiative cooling even in (rare) cases of nightime warming (e.g. warm front passing through) Modified: STULL -- Problem Area 4 - In case of several upper air stations, the upper air station selected for the computation of the lapse rate (IUPT) was not necessarily used for the computation of the lapse rate at the top of the boundary layer or the computation of the BG mixing height (the closest station was, whether IUPT or not) Modified: MIXHT, MIXHTST -- Problem Area 5 There is no check to make sure that the CALMET grid is located within the prognostic grid, thus possibly using far away and irrelevant prognostic gridpoints to initialize CALMET variables Modified: RDHD4, RDHD5, PARAMS.MET, READCF, INOUT (new) -- Problem Area 6 For NPSTA=-1 (precipitation interpolated from prognostic precipitation), the user-defined SIGMAP is not checked against the prgnostic gridsize thus allowing more than the 4 nearest progn. gridpoints to be used for the interpolatino onto the CALMET grid (possibly resulting in non-consistent precipitation fields) Modified: RDHD4, RDHD5, READHD ----------------------------------------------------------------------- ----------------------------------------------------------------------- 2. CALMET Changes (Version 5.815, Level 110421) ----------------------------------------------------------------------- *********************************************************************** --- Problem Area 1 --- When using multiple multi-hourly M3D files with overlapping times CALMET simulation stopped sooner than the requested end time. *********************************************************************** --------------------------------- Subroutine RDMM5 (Problem Area 1) --------------------------------- Ensure that the last record of the earlier M3D file is used and skip that hour in the next file. BEFORE change: -------------- if(mdathr.eq. (ndathr-isteppg+1) .and. ifirstpg.le.1) iskip=0 c --- Need something extra if run starts early on one day and the latest MM5 AFTER change: ------------- if(mdathr.eq. (ndathr-isteppg+1) .and. ifirstpg.le.1) iskip=0 c --- Make sure to read next record if overlapping files and multihourly c --- M3D timesteps (080710) if(mdathr.eq.ndathr.and. nfm3d.gt.1) iskip=1 c --- Need something extra if run starts early on one day and the latest MM5 *********************************************************************** --- Problem Area 2 --- When using convective overwater mixing height options (IMIXH>0) the overwater convective mixing height grows from the previous hour mixing height, not from the previous hour convective mixing height. *********************************************************************** --------------------------------- Subroutine WATER (Problem Area 2) --------------------------------- Save Convective mixing height array and pass on old values to MIXHMC and MIXHBG (rather than old values of the mixing height be it mechanical or convective) BEFORE change: -------------- real dcoast(mxnx,mxny),fcori(mxnx,mxny) AFTER change: ------------- real dcoast(mxnx,mxny),fcori(mxnx,mxny) real zicold(mxnx,mxny) BEFORE change: -------------- common /tjump/ dptt(mxnx,mxny) save zicold AFTER change: ------------- common /tjump/ dptt(mxnx,mxny) save zicold BEFORE change: -------------- zi(i,j)=amax1(zi(i,j),ziminw) if(imixh.eq.2) then c Batchvarova-Gryning convective mixing height call MIXHBG(ihrgmt,I,J,WT,DTHDZ,Tairow(nearow),Ilapse, : THRESHW,ZIMAXW,ZIMINW,USTAR(i,j),EL(i,j),ZI(i,j),ZICONV, : THT,THTP) else if (imixh.eq.1) then c Carson convective mixing height (similar to MIXHTST/MIXHSTS2) call MIXHMC(ihrgmt,i,j,wt,DTHDZ,ilapse,threshw, : zimaxw,zi(i,j),ziconv,dptt,THT,THTP) else c OCD mixing height (mechanical only) ziconv=0. endif AFTER change: ------------- zi(i,j)=amax1(zi(i,j),ziminw) zicold(i,j)=amin1(zicold(i,j),zi(i,j)) zicold(i,j)=amax1(zicold(i,j),ziminw) if(imixh.eq.2) then c Batchvarova-Gryning convective mixing height call MIXHBG(ihrgmt,I,J,WT,DTHDZ,Tairow(nearow),Ilapse, : THRESHW,ZIMAXW,ZIMINW,USTAR(i,j),EL(i,j),ZICOLD(i,j),ZICONV, : THT,THTP) Zicold(i,j)=ziconv else if (imixh.eq.1) then c Carson convective mixing height (similar to MIXHTST/MIXHSTS2) call MIXHMC(ihrgmt,i,j,wt,DTHDZ,ilapse,threshw, : zimaxw,zicold(i,j),ziconv,dptt,THT,THTP) Zicold(i,j)=ziconv else c OCD mixing height (mechanical only) ziconv=0. Zicold(i,j)=ziconv endif --------------------------------- Subroutine WATERP (Problem Area 2) --------------------------------- Save Convective mixing height array and pass on old values to MIXHMC and MIXHBG (rather than old values of the mixing height be it mechanical or convective) BEFORE change: -------------- real dcoast(mxnx,mxny),fcori(mxnx,mxny) AFTER change: ------------- real dcoast(mxnx,mxny),fcori(mxnx,mxny) real zicold(mxnx,mxny) BEFORE change: -------------- common /tjump/ dptt(mxnx,mxny) save zicold AFTER change: ------------- common /tjump/ dptt(mxnx,mxny) save zicold BEFORE change: -------------- zi(i,j)=amax1(zi(i,j),ziminw) if(imixh.eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG(ihrgmt,I,J,WT,DTHDZ,Tairow(nearow),1, : THRESHW,ZIMAXW,ZIMINW,USTAR(i,j),EL(i,j),ZI(i,j),ZICONV, : THT,THTP) else if (imixh.eq.1) then c --- Carson convective mixing height (similar to MIXHTST/MIXHSTS2) call MIXHMC(ihrgmt,i,j,wt,DTHDZ,1,threshw, : zimaxw,zi(i,j),ziconv,dptt,THT,THTP) else c OCD mixing height (mechanical only) ziconv=0. endif AFTER change: ------------- zi(i,j)=amax1(zi(i,j),ziminw) zicold(i,j)=amin1(zicold(i,j),zi(i,j)) zicold(i,j)=amax1(zicold(i,j),ziminw) if(imixh.eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG(ihrgmt,I,J,WT,DTHDZ,Tairow(nearow),1, : THRESHW,ZIMAXW,ZIMINW,USTAR(i,j),EL(i,j),ZICOLD(i,j),ZICONV, : THT,THTP) Zicold(i,j)=ziconv else if (imixh.eq.1) then c --- Carson convective mixing height (similar to MIXHTST/MIXHSTS2) call MIXHMC(ihrgmt,i,j,wt,DTHDZ,1,threshw, : zimaxw,zicold(i,j),ziconv,dptt,THT,THTP) Zicold(i,j)=ziconv else c OCD mixing height (mechanical only) ziconv=0. Zicold(i,j)=ziconv endif *********************************************************************** --- Problem Area 3 --- Radiative cooling always assumed at night when ITPROG>0 and lowest CALMET level is below lowest 3D.DAT level *********************************************************************** --------------------------------- Subroutine STULL (Problem Area 3) --------------------------------- When 3D.DAT lowest level is above lowest CALMET level(s) and ITPROG>0 allow the same nightime warming at lower CALMET levels as warming at the lowest 3D.DAT level (if any) ----------------------------------------------------------------------- BEFORE change: -------------- dt1= dts*exp(-alpha*z1/H) c assuming t1=t2 at sunset, t10=t20 => t1=dt1+t10=dt1+t20 t1= dt1+t20 AFTER change: ------------- dt1= dts*exp(-alpha*z1/H) c allow CALMET surface temp to warm up if prognostic near surface temp warm up if (dt2.gt.0) dt1=dt2 c assuming t1=t2 at sunset, t10=t20 => t1=dt1+t10=dt1+t20 t1= dt1+t20 *********************************************************************** --- Problem Area 4 --- Nearest Sounding was used for the computation of the lapse rate at the top of the boundary layer and for computation of BG mixing height rather than sounding IUPT (which may or may not be the same) *********************************************************************** --------------------------------- Subroutine MIXHT (Problem Area 4) --------------------------------- Use IUPT rather than IUSTA in all calls to MIXDT and MIXHBG ----------------------------------------------------------------------- BEFORE change: -------------- c --- Use 12z sounding to determine potential temperature lapse rate c --- in 'dzzi'-meter layer above last hour's convective mixing ht. htold=ziconv(i,j) nlev=nlbb(iusta) call mixdt(ihrgmt,i,j,iusta,htold,nlev,zlbb,tzbb,dptmin,dzzi, 1 tht,thtp,dtheta) if(abs(imixh).eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG (ihrgmt,I,J,WT,DTHETA,tzbb(iusta,1),2,THRESHL, : ZIMAX,ZIMIN,USTAR(i,j),EL(i,j),HTOLD,ZICONV(I,J), : THT,THTP) AFTER change: ------------- c --- Use 12z sounding to determine potential temperature lapse rate c --- in 'dzzi'-meter layer above last hour's convective mixing ht. htold=ziconv(i,j) nlev=nlbb(iupt) call mixdt(ihrgmt,i,j,iupt,htold,nlev,zlbb,tzbb,dptmin,dzzi, 1 tht,thtp,dtheta) if(abs(imixh).eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG (ihrgmt,I,J,WT,DTHETA,tzbb(iupt,1),2,THRESHL, : ZIMAX,ZIMIN,USTAR(i,j),EL(i,j),HTOLD,ZICONV(I,J), : THT,THTP) BEFORE change: -------------- nlev=nlaa(iusta) call mixdt(ihrgmt,i,j,iusta,htold,nlev,zlaa,tzaa,dptmin,dzzi, 1 tht,thtp,dtheta) if(abs(imixh).eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG (ihrgmt,I,J,WT,DTHETA,tzaa(iusta,1),2,THRESHL, : ZIMAX,ZIMIN,USTAR(i,j),EL(i,j),HTOLD,ZICONV(I,J), : tht,thtp) AFTER change: ------------- nlev=nlaa(iupt) call mixdt(ihrgmt,i,j,iupt,htold,nlev,zlaa,tzaa,dptmin,dzzi, 1 tht,thtp,dtheta) if(abs(imixh).eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG (ihrgmt,I,J,WT,DTHETA,tzaa(iupt,1),2,THRESHL, : ZIMAX,ZIMIN,USTAR(i,j),EL(i,j),HTOLD,ZICONV(I,J), : tht,thtp) --------------------------------- Subroutine MIXHTST (Problem Area 4) --------------------------------- Use IUPT rather than IUSTA in all calls to MIXDT and MIXHBG BEFORE change: -------------- c --- Use 12z sounding to determine lapse rate in 'dzzi'-meter layer c --- above last hour's convective mixing ht. htold=ziconv nlev=nlbb(iusta) call mixdt(ihrgmt,i,j,iusta,htold,nlev,zlbb,tzbb,dptmin,dzzi, 1 tht,thtp,dtheta) if(abs(imixh).eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG (ihrgmt,I,J,WT,DTHETA,tzbb(iusta,1),2,THRESHL, : ZIMAX,ZIMIN,USTAR,EL,HTOLD,ZICONV, : THT,THTP) AFTER change: ------------- c --- Use 12z sounding to determine lapse rate in 'dzzi'-meter layer c --- above last hour's convective mixing ht. htold=ziconv nlev=nlbb(iupt) call mixdt(ihrgmt,i,j,iupt,htold,nlev,zlbb,tzbb,dptmin,dzzi, 1 tht,thtp,dtheta) if(abs(imixh).eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG (ihrgmt,I,J,WT,DTHETA,tzbb(iupt,1),2,THRESHL, : ZIMAX,ZIMIN,USTAR,EL,HTOLD,ZICONV, : THT,THTP) BEFORE change: -------------- nlev=nlaa(iusta) call mixdt(ihrgmt,i,j,iusta,htold,nlev,zlaa,tzaa,dptmin,dzzi, 1 tht,thtp,dtheta) c if(abs(imixh).eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG (ihrgmt,I,J,WT,DTHETA,tzaa(iusta,1),2,THRESHL, : ZIMAX,ZIMIN,USTAR,EL,HTOLD,ZICONV, THT,THTP) AFTER change: ------------- nlev=nlaa(iupt) call mixdt(ihrgmt,i,j,iupt,htold,nlev,zlaa,tzaa,dptmin,dzzi, 1 tht,thtp,dtheta) c if(abs(imixh).eq.2) then c --- Batchvarova-Gryning convective mixing height call MIXHBG (ihrgmt,I,J,WT,DTHETA,tzaa(iupt,1),2,THRESHL, : ZIMAX,ZIMIN,USTAR,EL,HTOLD,ZICONV, THT,THTP) *********************************************************************** --- Problem Area 5 --- There is no check to make sure that the CALMET grid is located within the prognostic grid, thus possibly using far away and irrelevant prognostic gridpoints to initialize CALMET variables *********************************************************************** --------------------------------- Subroutine RDHD4 (Problem Area 5) --------------------------------- Check that MM4 gridpoints surround all CALMET gridpoints and print out QA3D.DAT file prior to that check BEFORE change: -------------- dimension neari(4),nearj(4),dnear(4) AFTER change: -------------- c --- CALMET corner points (110421) real*8 xc(4),yc(4) c --- MM4 corner points (defining MM4 boundary)(110421) real*8 xb(mxb3d),yb(mxb3d) dimension neari(4),nearj(4),dnear(4) BEFORE change: -------------- c --------------------------------------------------------------- c --- Find the 4 closest MM4 grid points to each CALMET grid point c (assume CALMET domain is inside MM4 grid section) c ---------------------------------------------------------------- AFTER change: -------------- c ----------------------------------------------------- c --- Print the MM4 grid points to the QA file c ----------------------------------------------------- if(lprt)then open(io4,file='QA3D.DAT',status='unknown') c write(io4,*)' MM4.DAT Grid Points' write(io4,*)' X Y Longitude Latitude' write(io4,*)' (km) (km) (deg E) (deg N)' do j=1,nyp do i=1,nxp write(io4,'(4f12.3)') xlcmm4(i,j),ylcmm4(i,j), 1 xlong4(i,j),xlat4(i,j) enddo enddo c close(io4) endif c ----------------------------------------------------------------------- c --- Ensure that the CALMET domain is within the MM4 sub-domain (110421) c ----------------------------------------------------------------------- c --- MM4 boundary points c --- Must include all boundary points and not just corner points c --- because of possible curvature of MM4 domain in CALMET coord. system do nc=1,nyp xb(nc)=xlcmm4(1,nc) yb(nc)=ylcmm4(1,nc) end do do i=1,nxp-1 nc=nyp+i xb(nc)=xlcmm4(1+i,nyp) yb(nc)=ylcmm4(1+i,nyp) end do do i=1,nyp-1 nc=nyp+nxp-1+i xb(nc)=xlcmm4(nxp,nyp-i) yb(nc)=ylcmm4(nxp,nyp-i) end do do i=1,nxp-1 nc=nyp+nxp+nyp-2+i xb(nc)=xlcmm4(nxp-i,1) yb(nc)=ylcmm4(nxp-i,1) end do c --- total number of boundary points nptbound=nyp+nxp+nyp+nxp-3 c --- CALMET corner points delg= dgrid * .001 xc(1) = xmap0 + (0.5 * delg) yc(1) = ymap0 + (0.5 * delg) xc(2) = xc(1) yc(2) = yc(1) + (ny-1) * delg xc(3) = xc(1) + (nx-1) * delg yc(3) = yc(1) xc(4) = xc(3) yc(4) = yc(2) c do nc=1,4 call inout(xc(nc),yc(nc),iflag,xb,yb,nptbound) if (iflag.eq.0) then c --- calmet gridpoint is outside of MM4 subdomain ->stop write(io6,*)'STOP - Subroutine RDHD4 -' write(io6,*)'CALMET grid outside of prognostic grid' write(6,*)'STOP - CALMET grid outside of prognostic grid' write(io6,*)'Check plot files QAMETG.BNA and QA3D.DAT' STOP end if end do c --------------------------------------------------------------- c --- Find the 4 closest MM4 grid points to each CALMET grid point c (assume CALMET domain is inside MM4 grid section) c ---------------------------------------------------------------- BEFORE change: -------------- c ----------------------------------------------------- c --- Print the MM4 grid points to the QA file c ----------------------------------------------------- if(lprt)then open(io4,file='QA3D.DAT',status='unknown') c write(io4,*)' MM4.DAT Grid Points' write(io4,*)' X Y Longitude Latitude' write(io4,*)' (km) (km) (deg E) (deg N)' do j=1,nyp do i=1,nxp write(io4,'(4f12.3)') xlcmm4(i,j),ylcmm4(i,j), 1 xlong4(i,j),xlat4(i,j) enddo enddo c close(io4) endif return end AFTER change: -------------- return end --------------------------------- Subroutine RDHD5 (Problem Area 5) --------------------------------- Check that MM5 gridpoints surround all CALMET gridpoints BEFORE change: -------------- data lprt/.true./ AFTER change: -------------- data lprt/.true./ c --- CALMET corner points (110421) real*8 xc(4),yc(4) c --- MM5 corner points (defining MM5 boundary)(110421) real*8 xb(mxb3d),yb(mxb3d) BEFORE change: -------------- c---------------------------------------------------------------- c --- Find the 4 closest MM5 grid points to each CALMET grid point c (assume CALMET domain is inside MM5 grid section) c --------------------------------------------------------------- AFTER change: -------------- c ---------------------------------------------------------------------- c --- Ensure that the CALMET domain is within the M3D sub-domain (110421) c ---------------------------------------------------------------------- c --- M3D boundary points: c --- Must include all boundary points and not just corner points c --- because of possible curvature of M3D domain in CALMET coord. system do nc=1,nyp xb(nc)=xlcmm4(1,nc) yb(nc)=ylcmm4(1,nc) end do do i=1,nxp-1 nc=nyp+i xb(nc)=xlcmm4(1+i,nyp) yb(nc)=ylcmm4(1+i,nyp) end do do i=1,nyp-1 nc=nyp+nxp-1+i xb(nc)=xlcmm4(nxp,nyp-i) yb(nc)=ylcmm4(nxp,nyp-i) end do do i=1,nxp-1 nc=nyp+nxp+nyp-2+i xb(nc)=xlcmm4(nxp-i,1) yb(nc)=ylcmm4(nxp-i,1) end do c --- total number of boundary points nptbound=nyp+nxp+nyp+nxp-3 c --- CALMET corner points delg= dgrid * .001 xc(1) = xmap0 + (0.5 * delg) yc(1) = ymap0 + (0.5 * delg) xc(2) = xc(1) yc(2) = yc(1) + (ny-1) * delg xc(3) = xc(1) + (nx-1) * delg yc(3) = yc(1) xc(4) = xc(3) yc(4) = yc(2) c do nc=1,4 call inout(xc(nc),yc(nc),iflag,xb,yb,nptbound) if (iflag.eq.0) then c --- calmet gridpoint is outside of M3D subdomain ->stop write(io6,*)'STOP - Subroutine RDHD5 -' write(io6,*)'CALMET grid outside of prognostic grid' write(6,*)'STOP - CALMET grid outside of prognostic grid' write(io6,*)'Check plot files QAMETG.BNA and QA3D.DAT' STOP end if end do c---------------------------------------------------------------- c --- Find the 4 closest MM5 grid points to each CALMET grid point c (assume CALMET domain is inside MM5 grid section) c --------------------------------------------------------------- --------------------------------- PARAMS.MET (Problem Area 5) --------------------------------- Define new parameter mxb3d (max number of gridpoints in MM5 subdomain) BEFORE change: -------------- parameter(mxwk3=mxwnd+2*mxnz+3) AFTER change: -------------- parameter(mxwk3=mxwnd+2*mxnz+3) parameter(mxb3d=2*mxnxp+2*mxnyp-3) --------------------------------- Subroutine READCF (Problem Area 5) --------------------------------- Create 'QAMETG.BNA' file outlining the CALMET boundaries to facilitate visual QA check (in conjunction with QA3D.DAT, containing the prognostic gridpoint locations in the same coordinate system) BEFORE change: -------------- c ----------------------------------------------------------- c ----------------- Print input parameters ------------------ c ----------------------------------------------------------- AFTER change: -------------- c ------------------------------------------------ c --- Print out Plot files with Met grid boundary c --- [BNA format] (110421) c ------------------------------------------------ open(io4,file='QAMETG.bna',status='REPLACE') c --- Set Upper Right met grid corner (km) c --- (Lower Left is xmap0,ymap0) xUR=xmap0+FLOAT(nx)*dgridkm yUR=ymap0+FLOAT(ny)*dgridkm write(io4,'(a15)') '" Met","grid",5' write(io4,'(f12.4,a1,f12.4)') xmap0,',',ymap0 write(io4,'(f12.4,a1,f12.4)') xmap0,',',yUR write(io4,'(f12.4,a1,f12.4)') xUR, ',',yUR write(io4,'(f12.4,a1,f12.4)') xUR, ',',ymap0 write(io4,'(f12.4,a1,f12.4)') xmap0,',',ymap0 close(io4) c ----------------------------------------------------------- c ----------------- Print input parameters ------------------ c ----------------------------------------------------------- --------------------------------- Subroutine INOUT (Problem Area 5) --------------------------------- Add new subroutine to check whether a gridpoint is within given boundaries at the end of the Fortran Code BEFORE change: -------------- 150 continue return end AFTER change: -------------- 150 continue return end c------------------------------------------------------------------------------ subroutine inout(xx,yy,iflag,xb,yb,nbpts) c------------------------------------------------------------------------------ c c --- CALMET Version: 5.815 Level: 110421 INOUT c --- Copied from V6.327 Level 090511 (F.Robe) c --- Based on BOUND_2.FOR, V2.0, J. Scire c c --- PURPOSE: performs in/out test on single point (xx,yy) for boundary c defined by xb,yb and returns in/out/on flag c c c --- UPDATES: c --- Bound_2.for v2.0 (2000-01-11) to V6.327 Level 090511 (F.Robe): c - Adapted to CALMET (variable definitions, conventions, c error messages,number of boundaries limited to 1 ) c c --- INPUT: c xx, YY - real*8 - coordinates of CALMET gridpoint (km) c XB (mxb3d),YB(mxb3d) - real*8 - Coordinates defining the boundary c the boundary must be closed (start c point=end point) c c via PARAMS.MET: io6,mxb3d (number of MM4/MM5 boundary points) c c NBPTS - integer - Number of points defining the boundary c --- OUTPUT: c IFLAG - integer - IN/OUT/ON flag c = 0: CALMET gridpoint OUTside of boundaries c = 1: CALMET gridpoint INside of boundaries c =-1: CALMET gridpoint ON boundaries c c --- INOUT called by: RDHD4, RDHD5 c --- INOUT calls: - c c ------------------------------------------------------------------------- include 'params.met' c --- CALMET gridpoint (dubbed "receptor" in the subroutine) real*8 xx,yy c --- boundary information real*8 xb(mxb3d),yb(mxb3d) real*8 small parameter(small=1.0d-3) !tolerance for being on boundary c --- test on boundary real*8 xs,ys,xp,yp,xe,ye,d2,ds,dc,small2 c --- work array with segment no. of all segments intersecting Y coordinate c --- of receptor integer ibsave(mxb3d) small2 = small*small c ----------------------------------------------- c --- STEP 1 - Find all boundary segments which c --- cross Y of the receptor c ----------------------------------------------- nseg=0 c --- Start and end indices of boundary points ns = 1 ne = nbpts c set nflag = 1 if we're starting above the line, =-1 if below c (do loop necessary if boundary is on y = yy) do j = ne,ns,-1 if (yb(j).gt.yy) then nflag = 1 goto 100 elseif (yb(j).lt.yy) then nflag = -1 goto 100 endif c if boundary is a straight line on yy, nflag = -1 nflag = -1 enddo 100 continue c for each boundary, make list of segments intersecting Y coordinate do j = ns+1,ne if(nflag.eq.1.and.yb(j).lt.yy) then nseg = nseg + 1 ibsave(nseg) = j nflag = -1 elseif(nflag.eq.-1.and.yb(j).gt.yy) then nseg = nseg + 1 ibsave(nseg) = j nflag = 1 endif c check if point is within 'small' of the boundary c boundary segment vector: xs = xb(j)-xb(j-1) ys = yb(j)-yb(j-1) ds = (xs*xs+ys*ys) !boundary segment length if (ds.ne.0.0) then c quietly ignore zero length segments c receptor vector: xp = xx-xb(j-1) yp = yy-yb(j-1) c perpendicular vector from boundary segment to receptor: xe = (xp - ((xs*xp+ys*yp)/ds)*xs) ye = (yp - ((xs*xp+ys*yp)/ds)*ys) d2 = (xe*xe+ye*ye) !perpendicular distance squared to receptor c projected position of perpendicular on segment: dc = (xs*xp+ys*yp) if (d2.lt.small2.and.dc.ge.(0.0-small2).and. & dc.le.(ds+small2)) then iflag = -1 goto 101 c return endif endif enddo c --- The number of segments intersecting y=yy must be zero or c --- an even number if(mod(nseg,2).ne.0)then write( 6,*)'ERROR -- STOP in subroutine INOUT' write(io6,*)'ERROR -- STOP in subroutine INOUT' write(io6,*)'Number of segments intersecting yy ', 1 'must be zero or an even number -- NSEG = ',nseg write(io6,*)'CALMET coordinates x,y:',xx,yy write(io6,*)'IBSAVE = ',(ibsave(n),n=1,nseg) stop endif c c --- If no boundary segments intersect yy, receptor must be OUTSIDE c --- of the boundary if(nseg.eq.0) then iflag = 0 goto 101 endif c c --------------------------------------------------- c --- STEP 2 - count no. of segments on either side of xx c --- even -> out c --- odd -> in c --- <> nseg -> boundary c --------------------------------------------------- c nlow = 0 nhigh = 0 do i = 1,nseg i1 = ibsave(i)-1 i2 = ibsave(i) xint = xb(i2)-(xb(i2)-xb(i1))*(yb(i2)-yy) & /(yb(i2)-yb(i1)) if (abs(xint-xx).lt.small) then c boundary: iflag = -1 goto 101 c return elseif (xint.gt.xx) then nhigh = nhigh + 1 elseif (xint.lt.xx) then nlow = nlow + 1 endif enddo if ((nhigh+nlow).eq.nseg) then if (mod(nhigh,2).eq.0) then c out: iflag = 0 else c in: iflag = 1 endif else c boundary: iflag = -1 endif 101 continue return end c --------------------------------------------------------------------- -- Problem Area 6 For NPSTA=-1 (precipitation interpolated from prognostic precipitation), the user-defined SIGMAP is not checked against the prgnostic gridsize thus allowing more than the 4 nearest progn. gridpoints to be used for the interpolatino onto the CALMET grid (possibly resulting in non-consistent precipitation fields) --------------------------------- Subroutine RDHD4 (Problem Area 6) --------------------------------- Compute MM4 gridsize and pass it through calling list BEFORE change: -------------- subroutine rdhd4 AFTER change: -------------- subroutine rdhd4(dxmm5) BEFORE change: -------------- 50 continue c AFTER change: -------------- 50 continue c c --- 110421- Compute MM5 gridsize in km dxmm5= sqrt ( (xlcmm4(1,1)-xlcmm4(2,1))**2 + : (ylcmm4(1,1)-ylcmm4(2,1))**2 ) c --------------------------------- Subroutine RDHD5 (Problem Area 6) --------------------------------- Compute MM5 gridsize and pass it via calling list BEFORE change: -------------- subroutine rdhd5(icloud,itwprog) AFTER change: -------------- subroutine rdhd5(icloud,itwprog,dxmm5) BEFORE change: -------------- return end AFTER change: -------------- c --- 110421- Compute MM5 gridsize in km dxmm5= sqrt ( (xlcmm4(1,1)-xlcmm4(2,1))**2 + : (ylcmm4(1,1)-ylcmm4(2,1))**2 ) return end --------------------------------- Subroutine READHD (Problem Area 6) --------------------------------- Check User-Selected SIGMAP against progn. gridsize and reset if necessary BEFORE change: -------------- call rdhd4 AFTER change: -------------- call rdhd4(dxmm5) BEFORE change: -------------- call rdhd5 (icloud,itwprog) AFTER change: -------------- call rdhd5 (icloud,itwprog,dxmm5) BEFORE change: -------------- return end AFTER change: -------------- c --- 110421- Check value of SIGMAP against DXMM5 when NPSTA=-1 c dxmm5 < sigmap < sqrt(2)* dxmm5 (allow 5% round-off) if (npsta.eq.-1) then dxmax= (2**0.5)*dxmm5 *1.05 if ((sigmap.lt.dxmm5).or.(sigmap.gt.dxmax)) then write(io6,*) 'WARNING ' write(io6,*) 'SIGMAP RESET to SQRT(2)* MM5gridsize ' sigmap=dxmax endif endif return end