program etatopo implicit none include "ecommons.spec" C integer IM,JM,IMJM,IMT,JMT,LM,LSM,JM2,IM1 integer IM,JM,LM,IMJM,LSM,IMT,JMT,IM1,JM2 include "parmeta" parameter (IMJM=IM*JM-JM/2,IMT=2*IM-1,JMT=JM/2+1) parameter (im1=im-1,jm2=jm-2) c integer*4 nsorc2,nexvlf,isrch,zeffl,zefflc,check1,check2,check3 integer*4 nssbsp,nssb,nssb2m,nssbp,nssbpp,kne,knw,kse,ksw integer*4 khl00,khh00,khl22,khh22,ninc real*4 dtr,pihf,a,cd,rvk2,rz0,dmin,hmin,phl,dlmdl,dphdl,rssb real*4 sumtxp,sumty,sumtyp,vk,wtamin,area,rarea,htdir,ndir real*4 sl,saob,shhob,hi,hob,wtarea,hghmn,hodir real*4 rtsubi,z0,hhdir,alah,sumtx,psbm,pav,psm,shs,hwelp Cmp real*4 wbd,sbd,tlm0d,ctph0,stph0,tph0d,dlmd,dphd integer ifill,jx,ix,mthdz0,imin,imax,is,jmin,jmax,llsb integer js,kssb,hmin1,hmax2,hmin3,nexvlk,nlob,ilob,nhghmn integer nlim,imn,nhmm,iex,ndirp,kndirt,kldirt,kwdirt,km,kndirs integer kldirs,kwdirs,ktdirs,nhused,m,ksdb,kspb,kssdbf,kssdb integer ndbmin,nsbmin,nsbhd,NINTC,nphd integer mssdbf,msdb,mspb,mssdb integer*4 LONNW,JTSUBH,NTSUBH,NTSUBI,ITSUBI, & NHREAD,issb,jssb,iqsb,jqsb,kqsb Cmp real*4 lata,latb,lona,lonb,tlm,tph,stph,ctph,sinphi, + glat(imjm),glon(imjm),fact,coslam,wb,sb,tph0,dlm + ,tdlm,dph Cmp real*4 tlat,tlon,r30sh,tslat,rtsubh,alm,aph,assb,r30s integer NDB(IMJM,4),npb(IMJM,4),NDPSB(4) real HSB(IMJM,4),PSB(IMJM,4) &, HS(4) &, AGTK(IMJM) real SM1(IMT,JMT),AGT(IMT,JMT) real*4 dlo,dla integer*4 ntiles,ISTART,IEND,LL C PARAMETER (NSSB=04) PARAMETER (NTSUBH=01) parameter (nsorc2=4) parameter (nssbsp=nssb*nssb+1) parameter (nssb2m=nssb*2-1) parameter (nssbp=nssb+1) parameter (nssbpp=nssbp+1) PARAMETER (PIHF=1.570796,DTR=0.01745329) parameter (ndbmin=1) parameter (nsbmin=1) parameter (dlo=1./120.,dla=1./120.) parameter (ntiles=18*36) C c arrays for subbox calculations real hssb(nssb,nssb,imjm) integer ndsb(nssb,nssb,imjm) real hssa(imjm) integer ndsa(imjm) c arrays for effective z0 real z0eff(imjm,4) real zha(imjm,4) real zas(imjm,4) integer isorc(imjm) real ao(4) real ha(4) real ds(4) real z0eff2(im,jm,4) real zha2(im,jm,4) real zas2(im,jm,4) integer msq(2,2) c arrays for statistics integer nsdb(102),nspb(102) integer nssdb(102) integer nssdbf(nssbsp) integer kndir(nssbsp,4,nsorc2) integer kldir(nssbsp,4,nsorc2) integer kwdir(nssbsp,4,nsorc2) integer kndirx(4),kldirx(4),kwdirx(4) real andirx(4),aldirx(4),awdirx(4) integer kisrch(3,4,nsorc2) integer klsrch(3,4,nsorc2) integer kwsrch(3,4,nsorc2) integer kountm(nsorc2) integer kountl(nsorc2) integer kountw(nsorc2) real zemean(4,nsorc2),zlmean(4,nsorc2),zwmean(4,nsorc2) real hamean(4,nsorc2),hlmean(4,nsorc2),hwmean(4,nsorc2) real asmean(4,nsorc2),almean(4,nsorc2),awmean(4,nsorc2) real zemin(4,nsorc2),zemax(4,nsorc2) real hamin(4,nsorc2),hamax(4,nsorc2) real asmin(4,nsorc2),asmax(4,nsorc2) integer izemin(4,nsorc2),izemax(4,nsorc2) integer ihamin(4,nsorc2),ihamax(4,nsorc2) integer iasmin(4,nsorc2),iasmax(4,nsorc2) integer jzemin(4,nsorc2),jzemax(4,nsorc2) integer jhamin(4,nsorc2),jhamax(4,nsorc2) integer jasmin(4,nsorc2),jasmax(4,nsorc2) integer kzemin(4,nsorc2),kzemax(4,nsorc2) integer khamin(4,nsorc2),khamax(4,nsorc2) integer khl2(jm),khh2(jm) integer kasmin(4,nsorc2),kasmax(4,nsorc2) c arrays for exteme values real exvl(nssbpp) real exvlc(nssb) c arrays for accumulating obstacle heights integer nadir(4) real htadir(4) real hhadir(4) real*4 hgtk(imjm),DUM2D(imt,jmt),TMP(IM,JM) . ,smk(imjm) . ,hgt(imt,jmt),sm(imt,jmt) . ,hbm2(imjm) . ,tfn(ntiles),tfx(ntiles),tln(ntiles),tlx(ntiles) . ,almx,almn,apmx,apmn . ,awlo,aelo,asla,anla . ,alt0d,almd,aphd,p,q . ,ald1,ald2,ald3,ald4,ald5,ald6,ald7,ald8 . ,apd1,apd2,apd3,apd4,apd5,apd6,apd7,apd8 c integer*2 iht(1200,1200) c integer*4 i,j,k,l,n,lat,lon,khl,khh . ,nlon,nlat,irecl . ,nwd,ned,nsd,nnd . ,nlo1,nlo2,nla1,nla2 . ,nrrd,ilast character*255 fname character*8 ctile(ntiles) character*3 alon character*2 alat character*1 lndir,ltdir,akm C NAMELIST /SEARESO/ seares data hwelp/99999./ DATA msq/4,1,2,3/ c logical last,inside c c *** Create topo for eta grid. c Uses 30 s global data obtained from ftp site edcftp.cr.usgs.gov c c *** Code obtained from U. of Athens and modified at FSL. c c_______________________________________________________________________________ c c *** Fill eta model common blocks. c call eta_commons write(6,*) 'back from eta_commons' , seares c *** Create eta topo file. c c=============================================================================== c c *** Create topo tile filenames and lat, lon extents. c At FSL topo tiles are 10x10 degrees with a file name indicating c the SW corner of the tile (e.g. U40N110W). c C Cnew (formerly parameter statements in oper) Cmp C rewind 15 C C read(15,SEARESO) CC C write(6,SEARESO) Cmp C ctph0=cosd(tph0d) C stph0=sind(tph0d) ctph0=cos(tph0d/57.2958) stph0=sin(tph0d/57.2958) KNE=IM KNW=IM-1 KSE=1-IM KSW=-IM KHL00=1 KHH00=IMJM KHL22=2*IM+1 KHH22=IMJM-2*IM NINC=2*IM-1 Cnew C n=0 do lat=-90,80,10 do lon=-180,170,10 n=n+1 if (lon .lt. 0) then lndir='W' else lndir='E' endif if (lat .lt. 0) then ltdir='S' else ltdir='N' endif write(alon,'(i3.3)') abs(lon) write(alat,'(i2.2)') abs(lat) ctile(n)='U'//alat//ltdir//alon//lndir tfn(n)=float(lat) tfx(n)=float(lat+10) tln(n)=float(lon) tlx(n)=float(lon+10) enddo enddo c c *** Find lat-lon extremes for eta grid. write(6,*) 'wbd, sbd ', wbd,sbd c C call rtll( wbd, sbd,tlm0d,ctph0,stph0,ald1,apd1) C call rtll( 0., sbd,tlm0d,ctph0,stph0,ald2,apd2) C call rtll(-wbd, sbd,tlm0d,ctph0,stph0,ald3,apd3) C call rtll( wbd, 0.,tlm0d,ctph0,stph0,ald4,apd4) C call rtll(-wbd, 0.,tlm0d,ctph0,stph0,ald5,apd5) C call rtll( wbd,-sbd,tlm0d,ctph0,stph0,ald6,apd6) C call rtll( 0.,-sbd,tlm0d,ctph0,stph0,ald7,apd7) C call rtll(-wbd,-sbd,tlm0d,ctph0,stph0,ald8,apd8) c C almx=max(ald1,ald2,ald3,ald4,ald5,ald6,ald7,ald8)+dlo c apmx=max(apd1,apd2,apd3,apd4,apd5,apd6,apd7,apd8)+dla C almn=min(ald1,ald2,ald3,ald4,ald5,ald6,ald7,ald8)-dlo C apmn=min(apd1,apd2,apd3,apd4,apd5,apd6,apd7,apd8)-dla call corners(im,jm,imjm,tph0d,tlm0d,dlmd,dphd,apmn,almn,apmx,almx) c print *,'ETA grid window:' print *,' Northwest =',apmx,almn print *,' Southeast =',apmn,almx print *,' ' C stop c c *** Boundary point heights are set to zero (where hbm2=0). c do j=1,jm khl2(j)=im*(j-1)-(j-1)/2+2 khh2(j)=im*j-j/2-1 enddo c do k=1,imjm hbm2(k)=0. enddo c do j=3,jm2 khl=khl2(j) khh=khh2(j) C write(6,*) 'khl, khh ', khl,khh do k=khl,khh hbm2(k)=1. enddo enddo c c *** Initialize height (hsb), ocean pts (psb), and total pts (ndb) within c each eta grid diamond to zero. c do m=1,4 do k=1,imjm hsb(k,m)=0. psb(k,m)=0. ndb(k,m)=0 enddo enddo c c *** Process tiles that cover eta grid domain. c do n=1,ntiles Cmp C IF (almn .le. tlx(n) .and. almx .ge. tln(n) .and. C . apmn .le. tfx(n) .and. apmx .ge. tfn(n)) then check1=0 check2=0 check3=0 if (apmn .le. tfx(n) .and. apmx .ge. tfn(n)) then IF (almn.le.tlx(n).and.almx.ge.tln(n)) check1=1 Cmp West extending into East if (almn .lt. -180 .and. tln(n).gt.0 .and. + (almn+360).le.tlx(n)) check2=1 Cmp East extending into West if (almx .gt. 180 .and. tlx(n) .lt. 0 .and. + almn.ge.tln(n) .and. (almx-360).ge.tln(n) ) check3=1 endif C write(6,*) 'lon bounds: this tile: ', almn,almx,tln(n),tlx(n) if (check1 .eq. 1 .or. check2 .eq. 1 .or. check3 .eq. 1) then Cmp write(6,*) 'inside for ', tlx(n),tfx(n),check1,check2,check3 c fname=topo_in(1:index(topo_in,' ')-1)//ctile(n) c c ********* Check that tile parameters match expected 10x10 degree size. c awlo=tln(n) aelo=tlx(n) asla=tfn(n) anla=tfx(n) nlon=int((aelo-awlo)/dlo)+1 nlat=int((anla-asla)/dla)+1 C Cmp conditional compilation for word length on DEC being different C CTO#ifdef DEC CTO irecl=nlon*0.5 CTO#else irecl=nlon*2 CTO#endif c C write(6,*) 'awlo,aelo,almn,almx: ', awlo,aelo,almn,almx Ctst nwd=(awlo-almn)/dlo if (almn .lt. -180 .and. awlo .gt. 0) then C write(6,*) 'here (1)' nwd=(awlo-(almn+360.))/dlo elseif (almn .gt. 0 .and. awlo .lt. 0) then nwd=(awlo-(almn-360.))/dlo else nwd=(awlo-almn)/dlo endif if (aelo .gt. 0 .and. almx .lt. 0) then C write(6,*) 'here (2)' ned=(aelo-(almx+360.))/dlo elseif (aelo .lt. 0 .and. almx .gt. 180) then C write(6,*) 'here (2a)' ned=(aelo-(almx-360.))/dlo C write(6,*) 'aelo, almx-360,ned ', aelo, almx-360, ned else ned=(aelo-almx)/dlo endif nsd=(asla-apmn)/dla nnd=(anla-apmx)/dla c if (nwd .le. 0) then nlo1=-nwd-4 else nlo1=1 endif c if (ned .ge. 0) then nlo2=nlon-ned+4 else nlo2=nlon endif c if (nnd .ge. 0) then nla1=nnd-4 else nla1=1 endif c if (nsd .le. 0) then nla2=nlat+nsd+4 else nla2=nlat endif c if (nlat .ne. 1200 .or. nlon .ne. 1200) then print *,'Unexpected tile size:',nlon,nlat print *,'Should be 1200 x 1200' print *,'Abort...' stop endif C print*,' ','(',nlo1,'-',nlo2,'),(',nla1,'-',nla2,')' c c ********* Open input topo file. c open(unit=1,file=fname,status='old' . ,form='unformatted',access='direct',recl=irecl/4 . ,err=900) l=index(fname//' ',' ')-1 Cmp print *,'reading topo data ',fname(1:l) print *,'reading topo data ',fname(l-8:l) Cmp write(6,*) ' ' c c ********* Process row by row. c Cno last=.false. Cno nrrd=nla1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Cmp alt0d=anla-dla/2.-(nrrd-1)*dla alt0d=tfx(n) assb=float(nssb) rssb=1.0/assb r30s=1.0/120.0 r30sh=r30s*0.5 RTSUBH=R30sh/(NTSUBh*1.0) Cnew do while (.not. last) do 710 nrrd=1,nlat read(1,rec=nrrd)(iht(i,nrrd),i=1,nlon) C if (nrrd .eq. nlat/2) write(6,*) (iht(I,NRRD),I=20,30) Cmp swap the topo data element by element #ifdef LITTLE do I=1,NLON call swapval(iht(I,nrrd)) enddo #endif 710 continue close(1) Cmp at this point the 1200 X 1200 block of data has been read do 7109 nrrd=1,nlat tlat=(ALT0D-(2*nrrd-1)*r30sh) if (tlat.gt.apmx) go to 7109 if (tlat.lt.apmn) go to 7109 do 7108 i=nlo1,nlo2 tlon=(awlo*1.0+(2*i-1)*(1./240.)) C if ((tlon.gt.blonem).and.(tlon.lt.blonw ).or. C 1 (tlon.gt.blone ).and.(tlon.lt.blonwp).or. C 2 (tlon.gt.blonep).and.(tlon.lt.blonwq)) go to 7108 Cnew DO 7107 jtsubh=1,ntsubh tslat=tlat+(jtsubh*2-ntsubh)*rtsubh ntsubi=ntsubh-ifix(ntsubh*(1.0-cos(dtr*tslat))) rtsubi=r30sh/(ntsubi*1.0) DO 7107 itsubi=1,ntsubi ALM=(tlon+(ITSUBi*2-ntsubi)*RTSUBi) APH=tslat nhread=nhread+1 ALM=ALM*DTR APH=APH*DTR call PQK_NEW(alm,aph,inside,p,q,k) IF (.NOT.INSIDE) GO TO 7107 c assign the subboxes based on the values of p and q issb = 1 + ifix(p*assb) jssb = 1 + ifix(q*assb) iqsb = 1 + ifix(p*2.0) jqsb = 1 + ifix(q*2.0) c kqsb=msq(iqsb,jqsb) if (NDB(k,kqsb).ge.0) then HSB(k,kqsb)=HSB(k,kqsb)+iht(i,nrrd) NDB(k,kqsb)=NDB(k,kqsb)+1 endif if (ndsb(issb,jssb,k).ge.0) then hssb(issb,jssb,k)=hssb(issb,jssb,k)+iht(i,nrrd) ndsb(issb,jssb,k)=ndsb(issb,jssb,k)+1 if ((issb.lt.1).or.(issb.gt.nssb).or. 1 (jssb.lt.1).or.(jssb.gt.nssb).or. 2 (k .lt.1).or.(k .gt.imjm)) then write(6,'('' data written outside the box'',4i6,2f7.3)') 1 issb,jssb,k ,k,p,q endif endif nhused=nhused+1 c 7107 continue 7108 continue 7109 continue Cmp endif on whole section? ENDIF Cmp enddo on number of tiles? enddo CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC write(6,*) 'done reading data?' write(6,*) 'nssbsp= ', nssbsp write(6,*) 'imjm= ', imjm c zero out counters for filling of grid boxes c do 4115 k=1,102 nsdb(k)=0 nspb(k)=0 nssdb(k)=0 4115 continue do 4117 k=1,nssbsp nssdbf(k)=0 4117 continue c convert positive numbers of points, fill counters c do 5177 k=1,imjm do 5163 m=1,4 if (ndb(k,m).gt.0) then if (mod(isorc(k),2).eq.0) then isorc(k)=isorc(k)+1 endif endif Ctest C if (ndb(k,m).ge.ndbmin) then C ndb(k,m)=-ndb(k,m) C endif Ctest ksdb=min(1+iabs(ndb(k,m)),102) kspb=min(1+iabs(npb(k,m)),102) nsdb(ksdb)=nsdb(ksdb)+1 nspb(kspb)=nspb(kspb)+1 5163 continue kssdbf=1 do 5167 j=1,nssb do 5165 i=1,nssb Ctest? C if (ndsb(i,j,k).ge.nsbmin) then C ndsb(i,j,k)=-ndsb(i,j,k) C endif Ctest? kssdb=min(1+iabs(ndsb(i,j,k)),102) nssdb(kssdb)=nssdb(kssdb)+1 if (ndsb(i,j,k).lt.0) kssdbf=kssdbf+1 5165 continue 5167 continue isorc(k)=isorc(k)*2 nssdbf(kssdbf)=nssdbf(kssdbf)+1 5177 continue c c report the numbers of boxes and subboxes with each number of points c write(6,'('' '')') c mssdbf=0 do 5403 k=1,nssbsp km=k-1 write(6,'('' boxes with'',i8,'' subboxes filled:'',i10)') 1 km,nssdbf(k) mssdbf=mssdbf+nssdbf(k) 5403 continue write (6,'('' '',8x,'' total boxes:'',i10)') 1 mssdbf c write(6,'('' '')') c msdb=0 mspb=0 mssdb=0 do 5405 k=1,102 km=k-1 if (km.eq.101) then akm='>' km=km-1 else akm=' ' endif write(6,'('' topo,land,subboxes with '',a1,i3, 1 '' points:'',3i8)') 2 akm,km,nsdb(k),nspb(k),nssdb(k) msdb=msdb+nsdb(k) mspb=mspb+nspb(k) mssdb=mssdb+nssdb(k) 5405 continue write (6,'('' topo,land,subboxes '',1x,3x, 1 '' totals:'',3i8)') 2 msdb,mspb,mssdb c write(6,'('' '')') C c c *** All rows of topo data have been processed. c Now calculate the heights (hgtk) and sea mask (smk). C dtr=acos(-1.)/180. tph0=tph0d*dtr wb=wbd*dtr sb=sbd*dtr dlm=dlmd*dtr dph=dphd*dtr tph=sb-dph tdlm=dlm+dlm C write(6,*) 'original values ' C write(6,*) tph0,wb,sb,dlm,dph,tph,dtr do j=1,jm khl=im*(j-1)-(j-1)/2+1 khh=im*j-j/2 c tlm=wb-tdlm+mod(j+1,2)*dlm tph=tph+dph stph=sin(tph) ctph=cos(tph) c do k=khl,khh tlm=tlm+tdlm sinphi=ctph0*stph+stph0*ctph*cos(tlm) glat(k)=asin(sinphi) coslam=ctph*cos(tlm)/(cos(glat(k))*ctph0) . -tan(glat(k))*tan(tph0) coslam=min(coslam,1.) fact=1. if (tlm .gt. 0.0) fact=-1. glon(k)=-tlm0d*dtr+fact*acos(coslam) glat(k)=glat(k)/dtr glon(k)=-glon(k)/dtr C write(6,*) glat(k),glon(k) enddo enddo CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c do 1077 M=1,4 do 1077 k=1,imjm if (NDB(K,M) .le. 0) then C write(6,*) 'K,M,NDB(K,M)= ', K,M,NDB(K,M) endif C IF (NDB(K,M).GT.0) GO TO 131 C NO DATA IN THE SUB-BOX: ASSUME PERCENTAGE OF WATER C SURFACE OF THE SUB-BOX TO BE 50 PERCENT, AND LEAVE C ELEVATION ZERO: if (npb(k,m).le.0) then PSB(K,M)=50. else PSB(K,M)=PSB(K,M)/npb(K,M) endif IF (NDB(K,M).LE.0) THEN HSB(K,M)=0.0 ELSE HSB(K,M)=HSB(K,M)/NDB(K,M) Cmp write(6,*) 'ave hgt= ', HSB(K,M),NDB(K,M) ENDIF 1077 continue c CNEW DO 108 K=1,IMJM 108 HGTK(K)=0. C DO 109 K=KHL22,KHH22 C PSM=0. DO 111 M=1,4 PSM=PSM+PSB(K,M) 111 CONTINUE PAV=PSM*0.25 IF (PAV.LE.50.) GO TO 121 C (IF NO DATA WAS AVAILABLE WITHIN ANY OF THE FOUR SUB-BOXES, C OR THE PERCENTAGE OF WATER SURFACE IS EXACTLY 50, THE POINT C IS DECLARED AS LAND) C C THE POINT IS DEFINED AS SEA/LAKE: C*** CALCULATE THE ELEVATION BY ASSIGNING TO THE FOUR SUB-BOX SQUARE C*** THE ELEVATION OF THE SUB-BOX WHICH HAD THE MAXIMUM PERCENTAGE OF C*** WATER COVERAGE, PSB(K,M). (NOTE: "LARGE LAKES OR INLAND SEAS C*** WILL NOT BE CODED AS 100" ON THE NAVY TAPE !!) C*** PSBM=50. DO 411 M=1,4 IF(NDB(K,M).EQ.0.OR.PSB(K,M).LE.PSBM) GO TO 411 HGTK(K)=HBM2(K)*HSB(K,M) PSBM = PSB(K,M) 411 CONTINUE IF (HGTK(K).LE.HWELP) GO TO 109 WRITE(6,412) K,PSBM,HGTK(K) 412 FORMAT(' K=',I5,' PSBM=',F5.1,' HGTK(K)=',F5.0) C WATER POINT IS HIGHER THAN HWELP. DECLARE IT LAND: SMK(K)=0. GO TO 109 C C THE POINT IS DEFINED AS LAND: C 121 SMK(K)=0. C write(6,*) 'ever here?????' C HS (1)=AMAX1(HSB(K,4),HSB(K,1)) NDPSB(1)= NDB(K,4)+NDB(K,1) HS (2)=AMAX1(HSB(K,3),HSB(K,2)) NDPSB(2)= NDB(K,3)+NDB(K,2) HS (3)=AMAX1(HSB(K,4),HSB(K,3)) NDPSB(3)= NDB(K,4)+NDB(K,3) HS (4)=AMAX1(HSB(K,1),HSB(K,2)) NDPSB(4)= NDB(K,1)+NDB(K,2) C C NOW CALCULATE THE NUMBER OF PAIRS OF SUB-BOXES WHICH HAVE C NAVY MOUNTAIN DATA IN AT LEAST ONE OF THE TWO SUB-BOXES C C HOW MANY PAIRS OF SUB-BOXES ARE THERE WITH HEIGHT DATA? NPHD=4 DO 151 N=1,4 IF(NDPSB(N).EQ.0) NPHD=NPHD-1 151 CONTINUE IF(NPHD.GT.0) GO TO 112 C C THERE WAS NO DATA WITHIN THE GRID BOX AND ITS ELEVATION IS C LEFT ZERO. DON'T WORRY ABOUT IT: IF IT'S A LAND POINT ITS C ELEVATION WILL BE RAISED IN 'PTETA' TO MATCH THAT UNDER C THE LOWEST NEIGHBORING NON-ZERO WIND POINT. HOWEVER: C SHOULD IT BE A WATER POINT AFTER ALL? DECLARE IT A WATER C POINT IF THE AVERAGE OF THE FOUR SUB-BOXES NEIGHBORING TO C K AND SOUTH OF IT IS MORE THAN 50 PERCENT WATER: PAV=0.25*(PSB(K+KSE,3)+PSB(K+KSE,2) & +PSB(K+KSW,1)+PSB(K+KSW,2)) Cmp IF(PAV.GT.50.) SMK(K)=1. GO TO 109 C 112 SHS=0. DO 110 N=1,4 SHS=SHS+HS(N) 110 CONTINUE HGTK(K)=HBM2(K)*SHS/NPHD C 109 CONTINUE C****SMK**smk************************ C C 2700,1350 for 8 minute C 5400,2700 for 4 minute C C Use the eta_sea_us call if domain is within box extending from C 10N,180W to 80N,30W. This uses 2 minute data C C if (seares .eq. 2) then call eta_sea_us(smk) elseif (seares .eq. 4) then CALL ETA_SEA(SMK,5400,2700) elseif (seares .eq. 8) then CALL ETA_SEA(SMK,2700,1350) else write(6,*) 'BAD VALUE FOR SEARES!!!!!! ' , seares STOP endif C C************************************* C C FOR AESTHETIC REASONS, MAKE SURE POINTS ALONG THE BUFFER LINES C REMAIN DECLARED WATER DO 431 K=KHL22,KHH22 SMK(K)=1.+HBM2(K)*(SMK(K)-1.) 431 CONTINUE call conh12t(smk,imjm,1,dum2d,imt,jmt) write(6,*) 'sea mask ' do J=JMT,1,-2 C write(6,273)(dum2d(I,J),I=1,IMT,4) enddo C----------------------------------------------------------- C*** CALCULATE AND SAVE AVERAGE ELEVATIONS C*** DO 301 K=KHL00,KHH00 AGTK(K)=0. 301 CONTINUE DO 302 K=KHL00,KHH00 NSBHD=4 DO 303 N=1,4 IF(NDB(K,N).EQ.0)NSBHD=NSBHD-1 303 CONTINUE C if (NSBHD .ne. 0) write(6,*) 'NSBHD= ', NSBHD IF(NSBHD.EQ.0)GO TO 302 AGTK(K)=(HSB(K,1)+HSB(K,2)+HSB(K,3)+HSB(K,4))/ 1 NSBHD C write(6,*) 'AGTK(K)= ', AGTK(K) 302 CONTINUE C*** C*** THE silhouette/mean MOUNTAINS: AT LAND POINTS WITH CONCAVE C*** ACTUAL (AVERAGE) TOPOGRAPHY, REPLACE THE C*** SILHOUETTE ELEVATION WITH THE AVERAGE ELEVATION C*** (22 FEB 94) DO 311 K=KHL22,KHH22 IF(SMK(K).GT.0.5)GO TO 311 ALAH=AGTK(K+1)+AGTK(K+NINC)+AGTK(K-1) 1 +AGTK(K-NINC)+2.*(AGTK(K+KNE) 2 +AGTK(K+KNW)+AGTK(K+KSW)+AGTK(K+KSE)) 3 -12.*AGTK(K) IF(ALAH.GE.0.)HGTK(K)=AGTK(K) 311 CONTINUE c c If K is a point seen as a "valley" when looking at three-point c averages in any one of the four directions, and it is at the c same time not higher than each of its four nearest neighbors, c and also not higher than each of its four second-nearest c neighbors, choose always the mean elevation (23 June 97) c c (Saddle points are expected to be frequently declared mean as c a result) c DO 316 K=KHL22,KHH22 IF(SMK(K).GT.0.5)GO TO 316 if (agtk(k).gt.agtk(k+kne) .and. agtk(k).gt.agtk(k+ knw) .and. 2 agtk(k).gt.agtk(k+ksw) .and. agtk(k).gt.agtk(k+ kse)) 3 go to 316 if (agtk(k).gt.agtk(k+ 1) .and. agtk(k).gt.agtk(k+ninc) .and. 2 agtk(k).gt.agtk(k- 1) .and. agtk(k).gt.agtk(k-ninc)) 3 go to 316 sumtx =agtk(k- 1)+agtk(k)+agtk(k+ 1) sumtxp=agtk(k+ ksw)+agtk(k)+agtk(k+ kne) sumty =agtk(k-ninc)+agtk(k)+agtk(k+ninc) sumtyp=agtk(k+ kse)+agtk(k)+agtk(k+ knw) if (sumtx .lt.agtk(k+ ksw)+agtk(k-ninc)+agtk(k+ kse) .and. 2 sumtx .lt.agtk(k+ knw)+agtk(k+ninc)+agtk(k+ kne)) 3 hgtk(k)=agtk(k) if (sumtxp.lt.agtk(k-ninc)+agtk(k+ kse)+agtk(k+ 1) .and. 2 sumtxp.lt.agtk(k- 1)+agtk(k+ knw)+agtk(k+ninc)) 3 hgtk(k)=agtk(k) if (sumty .lt.agtk(k+ kse)+agtk(k+ 1)+agtk(k+ kne) .and. 2 sumty .lt.agtk(k+ ksw)+agtk(k- 1)+agtk(k+ knw)) 3 hgtk(k)=agtk(k) if (sumtyp.lt.agtk(k+ 1)+agtk(k+ kne)+agtk(k+ninc) .and. 2 sumtyp.lt.agtk(k-ninc)+agtk(k+ ksw)+agtk(k- 1)) 3 hgtk(k)=agtk(k) 316 CONTINUE c c the silhouette/mean topography is done write(6,'('' silhouette/mean topography is done'')') C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCC Z0EFF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C c process the high-resolution arrays for roughness data c c zero out the arrays for the grid box averages c do 4203 k=1,imjm hssa(k)=0.0 ndsa(k)=0 4203 continue c c invert the negative numbers and find grid-box average height c do 4303 k=1,imjm do 4302 jssb=1,nssb do 4301 issb=1,nssb if (ndsb(issb,jssb,k).lt.0) then ndsb(issb,jssb,k)=-ndsb(issb,jssb,k) endif if (ndsb(issb,jssb,k).gt.0) then hssa(k)=hssa(k)+hssb(issb,jssb,k) ndsa(k)=ndsa(k)+ndsb(issb,jssb,k) hssb(issb,jssb,k)=hssb(issb,jssb,k)/ 1 ndsb(issb,jssb,k) C write(6,*) 'found an hssb value! ', C + issb,jssb,hssb(issb,jssb,k) endif 4301 continue 4302 continue 4303 continue c fill in the blanks with the average over the whole box c do 4343 k=1,imjm ifill=0 if (ndsa(k).eq.0) then write(6,'('' We have problems, no heights in k'',i8,2f6.1)') 1 k,glat(k),glon(k) ndsa(k)=1 hssa(k)=0 else hssa(k)=hssa(k)/ndsa(k) endif do 4342 jssb=1,nssb do 4341 issb=1,nssb if (ndsb(issb,jssb,k).eq.0.0) then hssb(issb,jssb,k)=hssa(k) ndsb(issb,jssb,k)=1 ifill=ifill+1 endif 4341 continue 4342 continue if (ifill.gt.0) then Cmp write(6,'('' '',i6,'' filled points for k='',i8)') Cmp 1 ifill,k endif 4343 continue c write(6,'('' empty subboxes have been filled'')') c do 4439 k=1,imjm C IMT=2*IM-1 JX=(K-1)/IMT+1 IX=K-(JX-1)*IMT IF(IX.LE.IM)THEN I=2*IX-1 J=2*JX-1 ELSE I=2*(IX-IM) J=2*JX ENDIF c write(6,'('' '',3i6,''=k,i,j'',f20.10,i10,''=h,n'')') c 1 k,i,j,hssa(k),ndsa(k) 4439 continue c c constants: c a=6371229.0 z0 = 0.15 c cd = 0.3 c cd = 0.8 cd = 0.6 vk = 0.4 rvk2 = 1.0 / ( vk ** 2 ) rz0 = 1.0 / z0 dmin=1.0e-6 hmin=2.0*z0+dmin wtamin=dmin c write(6,'('' constants defined'')') c c zero out statistics c do 4494 l=1,nsorc2 do 4493 m=1,4 hamean(m,l)=0.0 hlmean(m,l)=0.0 hwmean(m,l)=0.0 asmean(m,l)=0.0 almean(m,l)=0.0 awmean(m,l)=0.0 zemean(m,l)=0.0 zlmean(m,l)=0.0 zwmean(m,l)=0.0 hamin(m,l)= 1.0e30 hamax(m,l)=-1.0e30 asmin(m,l)= 1.0e30 asmax(m,l)=-1.0e30 zemin(m,l)= 1.0e30 zemax(m,l)=-1.0e30 ihamax(m,l)=0 ihamin(m,l)=0 iasmax(m,l)=0 iasmin(m,l)=0 izemax(m,l)=0 izemin(m,l)=0 jhamax(m,l)=0 jhamin(m,l)=0 jasmax(m,l)=0 jasmin(m,l)=0 jzemax(m,l)=0 jzemin(m,l)=0 khamax(m,l)=0 khamin(m,l)=0 kasmax(m,l)=0 kasmin(m,l)=0 kzemax(m,l)=0 kzemin(m,l)=0 do 4491 k=1,nssbsp kndir(k,m,l)=0 kldir(k,m,l)=0 kwdir(k,m,l)=0 4491 continue do 4492 k=1,3 kisrch(k,m,l)=0 kwsrch(k,m,l)=0 klsrch(k,m,l)=0 4492 continue 4493 continue 4494 continue c write(6,'('' statistics zeroed out'')') c c calculate 4 z0's for each grid box c do 4557 k=1,imjm C write(6,*) 'working k= ', K c l=isorc(k)+1 c c find the local mesh lengths c C---------------------------------------------------------- C*** C*** CONVERT FROM ONE-DIMENSIONAL K VALUES ON THE E-GRID C*** TO I,J VALUES ON THE FILLED E-GRID. C*** C---------------------------------------------------------- C C IMT=2*IM-1 JX=(K-1)/IMT+1 IX=K-(JX-1)*IMT IF(IX.LE.IM)THEN I=2*IX-1 J=2*JX-1 ELSE I=2*(IX-IM) J=2*JX ENDIF phl = ( j-1 ) * dphd + sbd c dlmdl = dlmd * dtr * cos(phl*dtr) dphdl = dphd * dtr c c find the length of the sides of the obstacles c ds(1)=a*rssb*dphdl ds(3)=a*rssb*dlmdl ds(2)=sqrt(ds(1)**2+ds(3)**2) ds(4)=ds(2) area=dphdl*dlmdl*a*a*2.0 rarea=1.0/area c if ((i.eq.j).or.((im-i).eq.(jm-j))) then Cmp write(6,'('' '',i6,''=i'',i6,''=j'')') i,j Cmp write(6,'('' ds ew nesw ns nwse'',4f12.3)') Cmp 1 ds endif C write(6,*) 'here 1' c c c test whether this is a boundary point c if (hbm2(k).eq.0.0) then c c set boundary points to z0 c do 4007 m=1,4 z0eff(k,m)=z0 zha(k,m)=0.0 zas(k,m)=0.0 4007 continue C write(6,*) 'setting boundary point' else c c calculate z0eff for interior points c mthdz0=2 c c mthdz0 1, find obstacles by searching for c relative minima and maxima along lines in each direction c Abandoned, as this was an older and an inferior option, having had c no consolidation of obstacles. FM, June 97 c c mthdz0 2, find obstacles by searching for c relative minima and maxima along lines in each direction c and eliminating insignificant relative extrema c if (mthdz0.eq.2) then c c loop for directions c do 4287 m=1,4 hhdir=0.0 htdir=0.0 ndir=0.0 imin=1 if (mod(m,2).eq.1) then c cardinal direction, diagonal on grid imax=nssb2m else c diagonal direction, along grid imax=nssb endif c c loop to select a line of subboxes c do 4285 is=imin,imax if (mod(m,2).eq.1) then c cardinal direction, diagonal on grid jmin=max(0,is-nssb) jmax=min(nssbp,is+1) else c diagonal direction, along grid jmin=0 jmax=nssbp endif llsb=jmax-jmin+1 nexvlf=0 c c loop to search for extrema along the line of subboxes c do 4223 js=jmin,jmax c select subbox indices by direction indicator if (m.eq.1) then c east to west issb=1+is-js jssb=js kssb=k endif if (m.eq.2) then c southwest to northeast issb=js jssb=is kssb=k endif if (m.eq.3) then c south to north issb=js jssb=nssb-is+js kssb=k endif if (m.eq.4) then c southeast to northwest issb=is jssb=js kssb=k endif c reindex to the interior of a box if (issb.lt.1) then issb=issb+nssb kssb=kssb+ksw endif if (jssb.lt.1) then jssb=jssb+nssb kssb=kssb+kse endif if (issb.gt.nssb) then issb=issb-nssb kssb=kssb+kne endif if (jssb.gt.nssb) then jssb=jssb-nssb kssb=kssb+knw endif c find minima and maxima if (js.eq.jmin) then c first point, taken to be first minimum hmin1=hssb(issb,jssb,kssb) isrch=1 else c any other point if (isrch.eq.1) then c looking for the first minimum if (hmin1.gt.hssb(issb,jssb,kssb)) then c lower the first minimum hmin1=hssb(issb,jssb,kssb) elseif (hmin1.lt.hssb(issb,jssb,kssb)) then c we have a minimum, now look for maximum nexvlf=nexvlf+1 exvl(nexvlf)=hmin1 hmax2=hssb(issb,jssb,kssb) isrch=2 endif elseif (isrch.eq.2) then c looking for a maximum if (hmax2.lt.hssb(issb,jssb,kssb)) then c raise the maximum hmax2=hssb(issb,jssb,kssb) elseif (hmax2.gt.hssb(issb,jssb,kssb)) then c we have a maximum, now look for minimum hmin3=hssb(issb,jssb,kssb) isrch=3 endif elseif (isrch.eq.3) then c looking for a minimum if (hmin3.gt.hssb(issb,jssb,kssb)) then c lower the minimum hmin3=hssb(issb,jssb,kssb) endif if ((hmin3.lt.hssb(issb,jssb,kssb)).or. . (js.eq.jmax)) then c we either have another minimum or are done c add an obstacle to the totals c hodir=hmax2-0.5*(hmin1+hmin3) c htdir=htdir+hodir c hhdir=hhdir+hodir**2 c ndir=ndir+1 c add an obstacle to the array of extrema nexvlf=nexvlf+1 exvl(nexvlf)=hmax2 nexvlf=nexvlf+1 exvl(nexvlf)=hmin3 c move the minimum value to front of next mt hmin1=hmin3 c now we are looking for another maximum hmax2=hssb(issb,jssb,kssb) isrch=2 endif endif endif 4223 continue kisrch(isrch,m,1)=kisrch(isrch,m,1)+1 kisrch(isrch,m,l)=kisrch(isrch,m,l)+1 if (smk(k).gt.0.5) then kwsrch(isrch,m,1)=kwsrch(isrch,m,1)+1 kwsrch(isrch,m,l)=kwsrch(isrch,m,l)+1 else klsrch(isrch,m,1)=klsrch(isrch,m,1)+1 klsrch(isrch,m,l)=klsrch(isrch,m,l)+1 endif c c consolidation of obstacles c Cmp write(6,*) 'how many obstacles? ', nexvlf if (nexvlf.ge.5) then c c calculate z0eff of the line just processed, c assuming that the width of the obstacles is unity c c sl = length of the line along which obstacles are calculated c nlob = n of the peak of the last obstacle along the line c if (m.eq.1) then sl=llsb*ds(3)*2.0 elseif (m.eq.3) then sl=llsb*ds(1)*2.0 else sl=llsb*ds(2) endif c nlob=nexvlf-1 saob=0.0 shhob=0.0 do 4235 ilob=2,nlob,2 hi=exvl (ilob)-0.5*(exvl (ilob-1)+exvl (ilob+1)) saob=saob+hi shhob=shhob+hi**2 4235 continue if (saob.gt.hmin) then hob=shhob/saob else hob=hmin saob=hmin shhob=hmin**2 endif c c calculate a tentative value of z0eff for this line c if (0.5*rz0*hob .lt. 0) write(6,*) 'neg alog', 0.5*rz0*hob wtarea= 2 0.5*cd*rvk2*saob/sl + 3 1.0/((alog(0.5*rz0*hob))**2) if (wtarea.lt.wtamin) then write(6,'('' nonpositive tentative wtarea'', 2 f20.10,2i8)') 3 wtarea,k,m wtarea=wtamin endif zeffl =0.5*hob/exp(1.0/sqrt(wtarea)) Cmp write(6,*) 'in z0eff code, wtarea,zeffl = ', wtarea,zeffl c c of the one or more inside minima, find the highest c 4243 continue hghmn=exvl(3) nhghmn=3 nlim=3 c c nlim = n of the last minimum c c nlim = 5 if there is only one inside minimum; it is the highest c if(nexvlf.gt.5) then nlim=nexvlf-2 do 4245 imn=5,nlim,2 if (hghmn.lt.exvl(imn)) then hghmn=exvl(imn) nhghmn=imn endif 4245 continue endif c c the highest inside minimum has been found. c does removing it and compressing exvl increase z0eff of the line? c nhmm=nhghmn-1 do 4251 iex=1,nhmm exvlc(iex)=exvl(iex) 4251 continue if (exvl(nhghmn+1).gt.exvl(nhmm)) then exvlc(nhmm)=exvl(nhghmn+1) endif c c the higher of the two maxima has been selected and stored c now transfer the remainder of exvl to exvlc c do 4253 iex=nhghmn,nlim exvlc(iex)=exvl(iex+2) 4253 continue c c calculate z0eff of the compressed array, zefflc c nlob=nlim-1 saob=0.0 shhob=0.0 do 4255 ilob=2,nlob,2 hi=exvlc(ilob)-0.5*(exvlc(ilob-1)+exvlc(ilob+1)) saob=saob+hi shhob=shhob+hi**2 4255 continue if (saob.gt.hmin) then hob=shhob/saob else hob=hmin saob=hmin shhob=hmin**2 endif c c calculate a tentative value of z0eff for this line c wtarea= 2 0.5*cd*rvk2*saob/sl + 3 1.0/((alog(0.5*rz0*hob))**2) if (wtarea.lt.wtamin) then write(6,'('' nonpositive tentative wtarea'', 2 f20.10,2i8)') 3 wtarea,k,m wtarea=wtamin endif zefflc=0.5*hob/exp(1.0/sqrt(wtarea)) c if (zeffl.lt.zefflc) then c c switch to compressed exvl and either exit the consolidation code c or search for another minimum to remove c do 4257 iex=nhmm,nlim exvl(iex)=exvlc(iex) 4257 continue nexvlf=nlim c if (nexvlf.lt.5) go to 4263 c zeffl =zefflc go to 4243 else go to 4263 endif c endif c 4263 continue c c There has been at most only one obstacle to start with, c or consolidation of obstacles is completed. c Add obstacles, if any, to the totals c if (nexvlf.ge.3) then nlob=nexvlf-1 do 4267 ilob=2,nlob,2 hodir=exvl(ilob)-0.5*(exvl(ilob-1)+exvl(ilob+1)) htdir=htdir+hodir hhdir=hhdir+hodir**2 ndir=ndir+1 4267 continue endif 4285 continue c in case we didn't find any obstacles if (ndir.eq.0) then ndir=1 htdir=hmin hhdir=hmin**2 endif c fill directional obstacle totals nadir(m)=ndir htadir(m)=htdir hhadir(m)=hhdir ndirp=ndir+1 kndir(ndirp,m,1)=kndir(ndirp,m,1)+1 kndir(ndirp,m,l)=kndir(ndirp,m,l)+1 if (smk(k).gt.0.5) then kwdir(ndirp,m,1)=kwdir(ndirp,m,1)+1 kwdir(ndirp,m,l)=kwdir(ndirp,m,l)+1 else kldir(ndirp,m,1)=kldir(ndirp,m,1)+1 kldir(ndirp,m,l)=kldir(ndirp,m,l)+1 endif 4287 continue c c end of mthdz0 2 c else c c default mthdz0, all h and a are based on hmin c do 4428 m=1,4 nadir(m)=1 htadir(m)=hmin hhadir(m)=hmin**2 4428 continue c c end of default mthdz0 c endif C write(6,*) 'here 2' c if ((i.eq.j).or.((im-i).eq.(jm-j))) then Cmp write(6,'('' na ew nesw ns nwse'',4i12)') C 1 nadir C write(6,'('' hta ew nesw ns nwse'',4f12.4)') C 1 htadir C write(6,'('' hha ew nesw ns nwse'',4f12.2)') Cmp 1 hhadir endif c do 4437 m=1,4 c c find the total area of obstacles in each direction c ao(m)=htadir(m)*ds(m)*rarea c c find the average heights of obstacles, use hmin to avoid zero c the average is weighted by the height of the obstacle c if (htadir(m) .lt. 0) write(6,*) 'divide by zero!' ha(m)=hhadir(m)/htadir(m) c 4437 continue C write(6,*) 'here 3' c c if ((i.eq.j).or.((im-i).eq.(jm-j))) then Cmp write(6,'('' ao ew nesw ns nwse'',4f12.5)') C 1 ao C write(6,'('' ha ew nesw ns nwse'',4f12.3)') Cmp 1 ha endif c c calculate z0eff c do 4552 m=1,4 zha(k,m)=ha(m) zas(k,m)=ao(m) c c accumulate statistics on zha and zas c if (hamin(m,1).gt.zha(k,m)) then hamin(m,1)=zha(k,m) ihamin(m,1)=i jhamin(m,1)=j khamin(m,1)=k endif if (hamin(m,l).gt.zha(k,m)) then hamin(m,l)=zha(k,m) ihamin(m,l)=i jhamin(m,l)=j khamin(m,l)=k endif if (hamax(m,1).lt.zha(k,m)) then hamax(m,1)=zha(k,m) ihamax(m,1)=i jhamax(m,1)=j khamax(m,1)=k endif if (hamax(m,l).lt.zha(k,m)) then hamax(m,l)=zha(k,m) ihamax(m,l)=i jhamax(m,l)=j khamax(m,l)=k endif hamean(m,1)=hamean(m,1)+zha(k,m) hamean(m,l)=hamean(m,l)+zha(k,m) if (smk(k).gt.0.5) then hwmean(m,1)=hwmean(m,1)+zha(k,m) hwmean(m,l)=hwmean(m,l)+zha(k,m) else hlmean(m,1)=hlmean(m,1)+zha(k,m) hlmean(m,l)=hlmean(m,l)+zha(k,m) endif c if (asmin(m,1).gt.zas(k,m)) then asmin(m,1)=zas(k,m) iasmin(m,1)=i jasmin(m,1)=j kasmin(m,1)=k endif if (asmin(m,l).gt.zas(k,m)) then asmin(m,l)=zas(k,m) iasmin(m,l)=i jasmin(m,l)=j kasmin(m,l)=k endif if (asmax(m,1).lt.zas(k,m)) then asmax(m,1)=zas(k,m) iasmax(m,1)=i jasmax(m,1)=j kasmax(m,1)=k endif if (asmax(m,l).lt.zas(k,m)) then asmax(m,l)=zas(k,m) iasmax(m,l)=i jasmax(m,l)=j kasmax(m,l)=k endif asmean(m,1)=asmean(m,1)+zas(k,m) asmean(m,l)=asmean(m,l)+zas(k,m) if (smk(k).gt.0.5) then awmean(m,1)=awmean(m,1)+zas(k,m) awmean(m,l)=awmean(m,l)+zas(k,m) else almean(m,1)=almean(m,1)+zas(k,m) almean(m,l)=almean(m,l)+zas(k,m) endif c c test for out-of-range values c if (zha(k,m).le.hmin) then c write(6,'('' nonpositive zha(k,m,1)'',f20.10,2i8)') c 1 zha(k,m),k,m zha(k,m)=hmin endif c if (zas(k,m).lt.0.0) then write(6,'('' negative zha(k,m,1)'',f20.10,2i8)') 1 zas(k,m),k,m zas(k,m)=0.0 endif c c calculate the effective z0 c from zas and zha c wtarea= 1 0.5*cd*rvk2*zas(k,m)+ 2 1.0/((alog(0.5*rz0*zha(k,m)))**2) if (wtarea.lt.wtamin) then write(6,'('' nonpositive wtarea'',f20.10,2i8)') 1 wtarea,k,m wtarea=wtamin endif z0eff(k,m)=0.5*zha(k,m)/exp(1.0/sqrt(wtarea)) C if (z0eff(k,m).ne. 0.15) then C write(6,*) 'zoeff= ', z0eff(k,m) C endif c c accumulate statistics c if (zemin(m,1).gt.z0eff(k,m)) then zemin(m,1)=z0eff(k,m) izemin(m,1)=i jzemin(m,1)=j kzemin(m,1)=k endif if (zemin(m,l).gt.z0eff(k,m)) then zemin(m,l)=z0eff(k,m) izemin(m,l)=i jzemin(m,l)=j kzemin(m,l)=k endif if (zemax(m,1).lt.z0eff(k,m)) then zemax(m,1)=z0eff(k,m) izemax(m,1)=i jzemax(m,1)=j kzemax(m,1)=k endif if (zemax(m,l).lt.z0eff(k,m)) then zemax(m,l)=z0eff(k,m) izemax(m,l)=i jzemax(m,l)=j kzemax(m,l)=k endif zemean(m,1)=zemean(m,1)+z0eff(k,m) zemean(m,l)=zemean(m,l)+z0eff(k,m) if (smk(k).gt.0.5) then zwmean(m,1)=zwmean(m,1)+z0eff(k,m) zwmean(m,l)=zwmean(m,l)+z0eff(k,m) else zlmean(m,1)=zlmean(m,1)+z0eff(k,m) zlmean(m,l)=zlmean(m,l)+z0eff(k,m) endif c 4552 continue c kountm(1)=kountm(1)+1 kountm(l)=kountm(l)+1 if (smk(k).gt.0.5) then kountw(1)=kountw(1)+1 kountw(l)=kountw(l)+1 else kountl(1)=kountl(1)+1 kountl(l)=kountl(l)+1 endif c endif 4557 continue c write(6,'('' z0eff calculation done'')') c c print statistics c do 4876 l=nsorc2,1,-1 c c report the numbers of boxes with each number of obstacles write(6,'('' '')') write(6,'('' statistics for topography source data'',i4)') . l write(6,'('' '')') c kndirt=0 kldirt=0 kwdirt=0 do 4561 m=1,4 kndirx(m)=0 kldirx(m)=0 kwdirx(m)=0 andirx(m)=0.0 aldirx(m)=0.0 awdirx(m)=0.0 4561 continue do 4563 k=1,nssbsp km=k-1 kndirs=kndir(k,1,l)+kndir(k,2,l)+kndir(k,3,l)+kndir(k,4,l) kldirs=kldir(k,1,l)+kldir(k,2,l)+kldir(k,3,l)+kldir(k,4,l) kwdirs=kwdir(k,1,l)+kwdir(k,2,l)+kwdir(k,3,l)+kwdir(k,4,l) ktdirs=kndirs+kldirs+kwdirs if (ktdirs.gt.0) then write(6,'('' '',i3,''o'',4(1x,i6,i5,i6))') 1 km,(kndir(k,m,l),kldir(k,m,l),kwdir(k,m,l),m=1,4) do 4562 m=1,4 kndirx(m)=kndirx(m)+kndir(k,m,l) kldirx(m)=kldirx(m)+kldir(k,m,l) kwdirx(m)=kwdirx(m)+kwdir(k,m,l) andirx(m)=andirx(m)+km*kndir(k,m,l) aldirx(m)=aldirx(m)+km*kldir(k,m,l) awdirx(m)=awdirx(m)+km*kwdir(k,m,l) 4562 continue kndirt=kndirt+kndirs kldirt=kldirt+kldirs kwdirt=kwdirt+kwdirs endif 4563 continue do 4564 m=1,4 if (kndirx(m).gt.0) then andirx(m)=andirx(m)/kndirx(m) endif if (kldirx(m).gt.0) then aldirx(m)=aldirx(m)/kldirx(m) endif if (kwdirx(m).gt.0) then awdirx(m)=awdirx(m)/kwdirx(m) endif 4564 continue do 4587 m=1,4 if (kountm(l).gt.0) then hamean(m,l)=hamean(m,l)/kountm(l) asmean(m,l)=asmean(m,l)/kountm(l) zemean(m,l)=zemean(m,l)/kountm(l) endif if (kountl(l).gt.0) then hlmean(m,l)=hlmean(m,l)/kountl(l) almean(m,l)=almean(m,l)/kountl(l) zlmean(m,l)=zlmean(m,l)/kountl(l) endif if (kountw(l).gt.0) then hwmean(m,l)=hwmean(m,l)/kountw(l) awmean(m,l)=awmean(m,l)/kountw(l) zwmean(m,l)=zwmean(m,l)/kountw(l) endif 4587 continue 4876 continue call conh12t(AGTK,imjm,1,AGT,imt,jmt) C write(6,*) 'zoeff values' do I=1,IMJM,50 C write(6,297) (z0eff(I,M),M=1,4) enddo 297 format(4(e10.4,x)) c do 8903 m=1,4 call conh12(z0eff(1,m),imjm,1,z0eff2(1,1,m),im,jm) 8903 continue write(6,*) '2d zoeff values' do J=JM,1,-2 C write(6,299) (z0eff2(I,J,2),I=1,IM,3) enddo 299 format(33(f4.2,x)) CCCCC write out the ZEFF file CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC open(unit=27,file='ZEFF',form='unformatted',access='sequential' + ,status='unknown') do m=1,4 tmp=-999. do j=1,jm do i=1,im TMP(I,J)=z0eff2(I,J,M) enddo enddo write(27) TMP enddo CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c c *** Reorder hgtk and smk array to old 2-d configuration. c call conh12t(hgtk,imjm,1,hgt,imt,jmt) call conh12t(smk,imjm,1,sm,imt,jmt) write(6,*) 'sea mask at end of topo ' do J=JMT,1,-JMT/10 write(6,273)(sm(I,J),I=1,IMT,IMT/10) enddo 273 format(80f2.0) C write(6,*) 'topo data before vertical interp (W,center,E) ' ISTART=1 do 276 LL=1,IMT/10 if (ISTART .ne. 1) ISTART=IEND+1 IEND=ISTART+10 C write(6,*) 'bounds are (LL,ISTART,IEND) ',LL, ISTART, IEND C write(6,*) ' ' istart=istart+1 276 continue do J=jmt,1,-4 write(6,275) (hgt(I,J),I=1,IMT,IMT/19) 275 format(22(f5.0,x)) enddo c *** Write eta grid topo. c open(1,file=topo_out,status='unknown',form='unformatted') write(1) hgt,sm close(1) stop c open(1,file='planine') c do j=jmt,1,-1 c write(1,'(1000F6.0)')(hgt(i,j),i=1,imt) c enddo c close(1) c Cstatic return c c *** Error trapping. c 900 print *,'Error opening input topo data.**********************' print *,' Filename = ',fname print*, '*************************************************' print*, '***********!!!!!!!!!!!!!!!!!********************' print*, '***********!!!!!!!!!!!!!!!!!********************' print*, '***********!!!!!!!!!!!!!!!!!********************' print*, ' M I S S I N G S O M E D A T A !!!!!!!!!!!!' print*, '***********!!!!!!!!!!!!!!!!!********************' print*, '***********!!!!!!!!!!!!!!!!!********************' print*, '***********!!!!!!!!!!!!!!!!!********************' print*, '*************************************************' stop c end c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& CC************************************************************ SUBROUTINE PQK_NEW &(ALM,APH,INSIDE,P,Q,K ) C ****************************************************************** C * * C * ROUTINE TO TRANSFORM ALAMBDA,APHI (LONGITUDE,LATITUDE) INTO * C * TRANSFORMED (ROTATED) LONGITUDE,LATITUDE, CHECK WHETHER THE * C * VALUES OBTAINED ARE WITHIN THE MODEL REGION, AND, IF THEY * C * ARE, CALCULATE VALUES OF P,Q (COORDINATES WITHIN A SQUARE * C * FORMED BY CONNECTING FOUR NEIGHBORING HEIGHT POINTS) AND K. * C * * C ****************************************************************** C * * C * GRID CONSTANTS: * C * WBD,SBD - WESTERN AND SOUTHERN BOUNDARIES ON LL OR TLL GRID * C * IN DEGREES * C * TLM0D,TPH0D - ANGLES OF ROTATION OF THE LL COORDINATE SYSTEM * C * IN THE DIRECTION OF LAMBDA AND PHI RESPECTIVELY * C * IN ORDER TO OBTAIN THE TLL COORDINATE SYSTEM * C * DLMD,DPHD - MESH SIDES IN DEGREES * C * * C ****************************************************************** c * 5/95 Wobus use whole domain, calculate the actual eta box k * c * and internal coordinates p and q here * c ****************************************************************** c include "parmeta" C include "parmgrd" include 'ecommons.h' C include "parmsub" C include "parmloc" PARAMETER (DTR=0.01745329,D50=0.50,H1=1.) common /pqkcnt/nwest,neast,nsouth,nnorth common /pqkcnt/kwest,keast,ksouth,knorth common /pqkcnt/kpneg,kqneg,kklow,kkhigh C---------------------------------------------------------------------- L O G I C A L & INSIDE C C---------------------------CONSTANTS----------------------------------- RIM1=IM1 TLM0=TLM0D*DTR TPH0=TPH0D*DTR DLM=DLMD*DTR DPH=DPHD*DTR RDLM=H1/DLM RDPH=H1/DPH crlwb c test for whole domain c ALMWB=WBD*DTR +DLM c APHSB=SBD*DTR +DPH c ALMEB=ALMWB +2*(IM-2)*DLM c APHNB=APHSB + (JM-3)*DPH c WB=ALMWB-DLM c SB=APHSB-DPH wb=wbd*dtr sb=sbd*dtr ALMWB=WB - DLM APHSB=SB - DPH ALMEB=ALMWB +(2*IM )*DLM APHNB=APHSB +( JM+1)*DPH crlwe STPH0= SIN(TPH0) CTPH0= COS(TPH0) C C--------------ALM,APH IS LON,LAT IN RADIAN MEASURE;-------------------- C------------TRANSFORM INTO ROTATED LONGITUDE,LATITUDE------------------ C RLM=ALM-TLM0 SRLM= SIN(RLM) CRLM= COS(RLM) SAPH= SIN(APH) CAPH= COS(APH) CC=CRLM*CAPH C ANUM=SRLM*CAPH DENOM=CTPH0*CC+STPH0*SAPH C ALM= ATAN2(ANUM,DENOM) APH= ASIN(CTPH0*SAPH-STPH0*CC) if (neast .lt. 5) then C write(6,*) 'bounds (w,s,e,n)' C write(6,*) almwb,aphsb,almeb,aphnb C write(6,*) 'we are at ',alm,aph endif C C-----IS ALM,APH INSIDE THE MODEL DOMAIN?------------------------------- C kmiss=0 if (alm.lt.almwb) then nwest=nwest+1 kmiss=kmiss-1 endif if (alm.gt.almeb) then neast=neast+1 kmiss=kmiss-2 endif if (aph.lt.aphsb) then nsouth=nsouth+1 kmiss=kmiss-4 endif if (aph.gt.aphnb) then nnorth=nnorth+1 kmiss=kmiss-8 endif C write(6,*) 'w,e,s,n ', nwest,neast,nsouth,nnorth if (kmiss.ne.0) go to 102 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc C C---------X1,Y1 IS A COORDINATE SYSTEM WITH DLM,DPH AS LENGTH UNITS----- c add y offset to index the centers of the boxes c instead of the southern corners X1=(ALM-WB)*RDLM Y1=(APH-SB)*RDPH+1.0 C---------X2,Y2 ROTATED FOR +45 DEG. & TRANSLATED FOR IM1--------------- X2=D50*( X1+Y1) Y2=D50*(-X1+Y1)+RIM1 C---------I2,J2 ARE COORDINATES OF center POINTS OF GRID BOXES-------- I2= INT(X2) J2= INT(Y2) C---------REMAINING PARAMETERS NEEDED TO KNOW THE POSITION-------------- C OF THE TRANSFORMED POINT WITHIN THE MODEL REGION P=X2-I2 Q=Y2-J2 C-----------------INDEX K CORRESPONDS TO I2,J2-------------------------- JR=J2-IM1 I3=I2-JR J3=I2+JR K=J3*IM-J3/2+(I3+2)/2 C---------------see if the point is in the domain----------------------- c note, i3 runs from 0 thru imt-1 c j3 runs from 0 thru jm-1 c P and q must be nonnegative for correct 0,0 point kmiss=0 if (i3.lt.0) then kwest=kwest+1 kmiss=kmiss+1 endif if (i3.ge.imt) then keast=keast+1 kmiss=kmiss+2 endif if (j3.lt.0) then ksouth=ksouth+1 kmiss=kmiss+4 endif if (j3.ge.jm) then knorth=knorth+1 kmiss=kmiss+8 endif if (p.lt.0) then kpneg=kpneg+1 kmiss=kmiss+16 endif if (q.lt.0.0) then kqneg=kqneg+1 kmiss=kmiss+32 endif if (k.le.0) then kklow=kklow+1 kmiss=kmiss+64 endif if (k.gt.imjm) then kkhigh=kkhigh+1 kmiss=kmiss+128 endif if (kmiss.gt.0) go to 102 c INSIDE=.TRUE. RETURN C----------------------------------------------------------------------- c 102 continue if (kmiss.eq.64) then write(6,'('' unexpected k too low '',3i8,2f10.6)') 1 k,i3,j3,p,q endif if (kmiss.eq.128) then write(6,'('' unexpected k too high'',3i8,2f10.6)') 1 k,i3,j3,p,q endif INSIDE=.FALSE. RETURN END c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine eta_sea(isea,idim,jdim) implicit none c include 'ecommons.h' c INTEGER IRES integer*4 nbx,nby,idim,jdim real*4 cxsta,cxsto,cysta,cysto,factor parameter (cxsta=-180.,cxsto=180.) parameter (cysta=-90.,cysto= 90.) c real*4 nbdnew,wbdnew,sbdnew,ebdnew . ,almin,almax,apmin,apmax . ,dummy,isea(IMJM) c integer*4 ista,isto,jsta,jsto . ,nx,ny . ,i,ii,j,jj,l,istatus c real seabyte(idim,jdim) nbx=idim nby=jdim write(6,*) 'starting sea program' c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 21600 minutes of longitude around earth IRES=21600/IDIM write(6,*) 'IRES= ', IRES factor=60./float(IRES) if (IRES .eq. 4) then write(6,*) 'opening 4 minute' open(unit=2,file='global_4m.ieee',form='unformatted', . access='sequential',err=920,status='old') elseif (IRES .eq. 8) then write(6,*) 'opening 8 minute' open(unit=2,file='global_8m.ieee',form='unformatted', . access='sequential',err=920,status='old') else write(6,*) 'opening something with an unsupported IRES value...' stop endif write(6,*) 'unit 2 open ', nbx,nby do J=1,nby C write(6,*) 'reading for J= ', J read(2) (seabyte(I,J),I=1,nbx) C write(6,*) 'read rec= ', J, (seabyte(I,J),I=20,1600,40) enddo goto 921 920 write(6,*) '!!!!!!!!!trouble opening the seamask data!!!!!!' 921 write(6,*) 'data read successfully ' close(2) c c *** Sector out a tile of sea from the global data set c that covers the ETA grid domain. write(6,*) 'wbd,sbd,tlm0d,tph0d= ',wbd,sbd,tlm0d,tph0d c C call rtll( 0,-sbd,tlm0d,ctph0,stph0,dummy,nbdnew) C call rtll( wbd,-sbd,tlm0d,ctph0,stph0,wbdnew,dummy) C call rtll(-wbd,-sbd,tlm0d,ctph0,stph0,ebdnew,dummy) C call rtll( wbd, sbd,tlm0d,ctph0,stph0,dummy,sbdnew) call corners(im,jm,imjm,tph0d,tlm0d,dlmd,dphd, + sbdnew,wbdnew,nbdnew,ebdnew) c almin=int(wbdnew/5)*5-5 almax=int(ebdnew/5)*5+5 apmin=int(sbdnew/5)*5-5 apmax=int(nbdnew/5)*5+5 c write(*,*)'nbdnew,wbdnew,ebdnew,sbdnew' write(*,*)nbdnew,wbdnew,ebdnew,sbdnew write(*,*)'apmax,almin,almax,apmin' write(*,*)apmax,almin,almax,apmin c write(6,*) 'cxsta= ', cxsta ista=int((almin-cxsta)*factor)+1 if (ista .lt. 0) then print*, ista, 'ista' C stop endif isto=int((almax-cxsta)*factor) C if (isto .lt. 0) isto=isto+nbx jsta=int((apmin-cysta)*factor)+1 jsto=int((apmax-cysta)*factor) nx=isto-ista+1 ny=jsto-jsta+1 c write(6,*) 'ista, isto, nx',ista,isto,nx write(6,*) 'jsta, jsto, ny',jsta,jsto,ny c call create_seatype(apmax,almin,ista,jsta + ,nx,ny,seabyte,nbx,nby,isea,IRES) return end c=============================================================================== c subroutine create_seatype(seanbd,seawbd,ista,jsta . ,irow,jcol,seabyte,nbx,nby,isea,IRES) c c *** Code to produce eta sea mask c Program written using the structure of mountnew.f c Written by F. Mesinger, W. Collins and Z. Janjic. c (Modified by S. Nickovic for surface types) c (Modified for use at FSL) c (modified for use at NCEP) c c *** Specify the following: c tlm0d,tph0d: transformed lambda,phi, in degrees c (origin of the rotated coordinate system) c wbd,sbd: western and southern boundary of the c transformed grid, in degrees of that grid c dlmd,dphd: mesh sides, in degrees c c dld - resolution of the input data c c *** Written by S. Nickovic, Univ. of Athens, Jan. 1996 c c Method: The largest percent of a seamask type that c enters in the h grid-box is used as c the representative one. c implicit none c include 'ecommons.h' c integer*4 ntex,irow,jcol,IRES real*4 dld parameter(ntex=2) c real*4 dum2d(imt,jmt) . ,alt0d,almd,aphd,p,q . ,seanbd,seawbd,lrat,tsum c integer*4 ndb(IMJM,ntex) . ,ilast,inimax,iindex . ,ista,jsta,nbx,nby . ,ksea,itmp,i,j,k,l,ii,jj,istatus c c logical last,inside c real*4 isea(imjm) real seabyte(nbx,nby),itypd(irow,jcol) c_______________________________________________________________________________ dld=float(IRES)/60. c c *** Read sea masks from ETA topo file. c c *** Initialize arrays. c do k=1,imjm isea(k)=99. enddo do i=1,ntex do k=1,imjm ndb(k,i)=0 enddo enddo c if (jsta .gt. nby) then write(6,*) 'big trouble!!!!' stop endif do j=jcol,1,-1 Csauditest do j=1,jcol jj=jsta+j-1 do i=1,irow ii=ista+i-1 if (ii .lt. 0) ii=nbx+ii-1 if (mod(II,400).eq.0 .and. mod(JJ,400) .eq. 0) then write(6,*) 'value at ', II,JJ, '-->', I,J endif itypd(i,j)=seabyte(ii,jj) enddo enddo Cmp write(6,*) 'what is itypd at this point', irow,jcol do J=jcol,1,-40 write(6,243) (itypd(I,J),I=1,irow,50) enddo 243 format(61f2.0) Cmp c c *** Here starts processing row by row. c last=.false. c write(6,*) 'seanbd,dld= ', seanbd, dld do while (.not. last) c C Cmp loop from north to south C do j=jcol,1,-1 Cmp alt0d=seanbd+(j-2)*dld+dld*0.5 alt0d=seanbd-(jcol-j-2)*dld+dld*0.5 C write(6,*) 'alt0d= ', alt0d ilast=0 c do i=1,irow c c ************ sea mask value at this point is itypd(i,j). c ksea=itypd(i,j)+1 if (ksea .ne. 1 .and. ksea .ne. 2) write(6,*) 'ksea= ', ksea c almd=seawbd+(i-1)*dld Cold aphd=seanbd-(j-1)*dld aphd=seanbd-(JCOL-j-1)*dld C if (mod(I,20).eq.0 .and. mod(J,20) .eq. 0) then C write(6,*) 'I,J,aphd,almd ', I,J,aphd,almd C endif call pqk(almd,aphd,inside,p,q,k) C if (mod(I,20).eq.0 .and. mod(J,20) .eq. 0) then C write(6,*) 'inside= ', inside C endif if (inside) then ilast=ilast+1 c c c h(k+imt) c c h(k+im1) h(k+im) c c h(k) c c if (p .gt. 0.5 .and. q .le. 0.5) then k=k+im elseif (p .gt. 0.5 .and. q .gt. 0.5) then k=k+imt elseif (p .le. 0.5 .and. q .gt. 0.5) then k=k+im1 endif ndb(k,ksea)=ndb(k,ksea)+1 endif c enddo last=(ilast .eq. 0 .and. alt0d .lt. tph0d) enddo last=.true. enddo c c *** All rows of land/sea data have been processed. c Now find land/sea which has max contribution in a c h-grid box c Cmp values will be ones or twos (1=land, 2=water) C Cmp for output (0=land, 1=water) do k=1,imjm tsum=ndb(k,1)+ndb(k,2) if (tsum .ne. 0) then lrat= float(ndb(k,2)) / tsum C write(6,*) k,lrat,ndb(k,1),ndb(k,2) if (lrat .gt. 0.5) then C write(6,*) 'generating land' isea(k)=0. else C write(6,*) 'generating sea' isea(k)=1. endif else write(6,*) 'no data ', k C isea=0. endif enddo Cmptest c c *** Fill in for any missing points. c (Do recursively until all missing values have been filled.) last=.false. do while (.not. last) last=.true. do k=1,imjm if (isea(k).lt.0. .or. isea(k).gt.1.) then last=.false. if (k .le. im) then !Bottom row if (k .eq. 1) then !Lower-left corner isea(k)=isea(k+im) elseif (k .eq. im) then !Lower-right corner isea(k)=isea(k+im1) else isea(k)=min(isea(k+im1),isea(k+im)) endif elseif (k .gt. imjm-im) then !Top row if (k .eq. imjm-im1) then !Upper-left corner isea(k)=isea(k-im1) elseif (k .eq. imjm) then !Upper-right corner isea(k)=isea(k-im) else isea(k)=min(isea(k-im1),isea(k-im)) endif else !Interior row if (mod(k,imt) .eq. 1) then !Left column, odd row isea(k)=min(isea(k-im1),isea(k+1),isea(k+im)) elseif (mod(k,imt) .eq. im) then !Right column, odd row isea(k)=min(isea(k-im),isea(k-1),isea(k+im1)) else isea(k)=min(isea(k+im1),isea(k+im) . ,isea(k-im1),isea(k-im)) endif endif endif enddo enddo Cmptest c c write(6,*) 'GOT TO THIS POINT' C open(1,file=sea_out,status='unknown',form='unformatted') C write(1) isea C close(1) call conh12t(isea,imjm,1,dum2d,imt,jmt) do J=JMT,1,-JMT/10 C write(6,273)(dum2d(I,J),I=1,IMT,IMT/10) enddo 273 format(71f2.0) c return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine eta_sea_us(isea) implicit none c include 'ecommons.h' c integer*4 nbx,nby real*4 cxsta,cxsto,cysta,cysto,factor parameter(factor=60./2.) parameter (cxsta=180.,cxsto=330.) parameter (cysta= 10.,cysto= 80.) C parameter ( nbx=int(abs(cxsto-cxsta)*factor)+1, C + nby=int(abs(cysta-cysto)*factor)+1 ) parameter(nbx=4501,nby=2101) c real*4 nbdnew,wbdnew,sbdnew,ebdnew . ,almin,almax,apmin,apmax . ,dummy,isea(IMJM) c integer*4 ista,isto,jsta,jsto . ,nx,ny . ,i,ii,j,jj,l,istatus c real seabyte(nbx,nby) write(6,*) 'starting sea program ', nbx,nby c open(2,file='US_2m_slm.ieee',form='unformatted', . access='sequential', . status='old',ERR=919) write(6,*) 'unit 2 open ' do J=1,nby read(2) (seabyte(I,J),I=1,nbx) enddo goto 901 919 write(6,*) 'trouble opening the sea mask file' 901 write(6,*) ' done reading data ' close(2) c c *** Sector out a tile of soil types from the global data set c that covers the ETA grid domain. write(6,*) 'wbd,sbd,tlm0d,tph0d= ',wbd,sbd,tlm0d,tph0d c call rtll( 0,-sbd,tlm0d,ctph0,stph0,dummy,nbdnew) call rtll( wbd,-sbd,tlm0d,ctph0,stph0,wbdnew,dummy) call rtll(-wbd,-sbd,tlm0d,ctph0,stph0,ebdnew,dummy) call rtll( wbd, sbd,tlm0d,ctph0,stph0,dummy,sbdnew) c almin=int(wbdnew/5)*5-5 almax=int(ebdnew/5)*5+5 if (almin .lt. 0.) almin=almin+360. if (almax .lt. 0.) almax=almax+360. apmin=int(sbdnew/5)*5-5 apmax=int(nbdnew/5)*5+5 c write(*,*)'nbdnew,wbdnew,ebdnew,sbdnew' write(*,*)nbdnew,wbdnew,ebdnew,sbdnew write(*,*)'apmax,almin,almax,apmin' write(*,*)apmax,almin,almax,apmin c ista=int((almin-cxsta)*factor)+1 C if (ista .lt. 0) ista=ista+nbx isto=int((almax-cxsta)*factor) C if (isto .lt. 0) isto=isto+nbx jsta=int((apmin-cysta)*factor)+1 jsto=int((apmax-cysta)*factor) nx=isto-ista+1 ny=jsto-jsta+1 c write(6,*) 'ista, isto, nx',ista,isto,nx write(6,*) 'jsta, jsto, ny',jsta,jsto,ny c call create_seaus(apmax,almin,ista,jsta + ,nx,ny,seabyte,nbx,nby,isea) return end c=============================================================================== c subroutine create_seaus(soilnbd,soilwbd,ista,jsta . ,irow,jcol,seabyte,nbx,nby,isea) c c *** Code to produce eta sea mask c Program written using the structure of mountnew.f c Written by F. Mesinger, W. Collins and Z. Janjic. c (Modified by S. Nickovic for surface types) c (Modified for use at FSL) c (modified for use at NCEP) c c *** Specify the following: c tlm0d,tph0d: transformed lambda,phi, in degrees c (origin of the rotated coordinate system) c wbd,sbd: western and southern boundary of the c transformed grid, in degrees of that grid c dlmd,dphd: mesh sides, in degrees c c dld - resolution of the input data c c *** Written by S. Nickovic, Univ. of Athens, Jan. 1996 c c Method: The largest percent of a seamask type that c enters in the h grid-box is used as c the representative one. c implicit none c include 'ecommons.h' c integer*4 ntex,irow,jcol real*4 dld parameter(ntex=2,dld=2./60.) c real*4 dum2d(imt,jmt) . ,alt0d,almd,aphd,p,q . ,soilnbd,soilwbd,lrat,tsum c integer*4 ndb(IMJM,ntex) . ,ilast,inimax,iindex . ,ista,jsta,nbx,nby . ,ksea,itmp,i,j,k,l,ii,jj,istatus c c logical last,inside c real*4 isea(imjm) real seabyte(nbx,nby),itypd(irow,jcol) c_______________________________________________________________________________ c c *** Read sea masks from ETA topo file. c write(6,*) 'inside subroutine...,imjm= ',imjm write(6,*) 'US data set ', nbx,nby c *** Initialize arrays. c do k=1,imjm isea(k)=-99. enddo do i=1,ntex do k=1,imjm ndb(k,i)=0 enddo enddo c if (jsta .gt. nby) then write(6,*) 'big trouble!!!!' stop endif Coirg do j=jcol,1,-1 do j=1,jcol jj=jsta+j-1 do i=1,irow ii=ista+i-1 if (mod(II,400).eq.0 .and. mod(JJ,400) .eq. 0) then write(6,*) 'value at ', II,JJ, '-->', I,J endif itypd(i,j)=seabyte(ii,jj) C if (seabyte(ii,jj) .lt. 0) itypd(i,j)=seabyte(ii,jj)+256 enddo enddo Cmp write(6,*) 'what is itypd at this point', irow,jcol do J=jcol,1,-40 write(6,243) (itypd(I,J),I=1,irow,50) enddo 243 format(61f2.0) Cmp c c *** Here starts processing row by row. c last=.false. c write(6,*) 'soilnbd,dld= ', soilnbd, dld do while (.not. last) c C Cmp loop from north to south C do j=jcol,1,-1 Cmp alt0d=soilnbd+(j-2)*dld+dld*0.5 alt0d=soilnbd-(jcol-j-2)*dld+dld*0.5 C write(6,*) 'alt0d= ', alt0d ilast=0 c do i=1,irow c c ************ sea mask value at this point is itypd(i,j). c ksea=itypd(i,j)+1 if (ksea .ne. 1 .and. ksea .ne. 2) write(6,*) 'ksea= ', ksea c almd=soilwbd+(i-1)*dld Cold aphd=soilnbd-(j-1)*dld aphd=soilnbd-(JCOL-j-1)*dld C if (mod(I,20).eq.0 .and. mod(J,20) .eq. 0) then C write(6,*) 'I,J,aphd,almd ', I,J,aphd,almd C endif call pqk(almd,aphd,inside,p,q,k) C if (mod(I,20).eq.0 .and. mod(J,20) .eq. 0) then C write(6,*) 'inside= ', inside C endif if (inside) then ilast=ilast+1 c c *************** The soiltexture (ksea) data point is within the box c formed by four neighboring h points with h(k) point c as south point. c c h(k+imt) c c h(k+im1) h(k+im) c c h(k) c c *************** Now calculate the number of points for each of 73 c soil types that are inside the h grid-box: c if (p .gt. 0.5 .and. q .le. 0.5) then k=k+im elseif (p .gt. 0.5 .and. q .gt. 0.5) then k=k+imt elseif (p .le. 0.5 .and. q .gt. 0.5) then k=k+im1 endif ndb(k,ksea)=ndb(k,ksea)+1 endif c enddo last=(ilast .eq. 0 .and. alt0d .lt. tph0d) enddo last=.true. enddo c c *** All rows of land/sea data have been processed. c Now find land/sea which has max contribution in a c h-grid box c Cmp values will be ones or twos (1=land, 2=water) C Cmp for output (0=land, 1=water) do k=1,imjm tsum=ndb(k,1)+ndb(k,2) if (tsum .ne. 0) then lrat= float(ndb(k,2)) / tsum C write(6,*) k,lrat,ndb(k,1),ndb(k,2) if (lrat .gt. 0.5) then C write(6,*) 'generating land' isea(k)=0. else C write(6,*) 'generating sea' isea(k)=1. endif else Cmp pure water write(6,*) 'no data' C isea=0. endif C write(6,*) ndb(k,1),ndb(k,2) enddo c write(6,*) 'past the do loop' Cmp recursively fill the missing values last=.false. do while (.not. last) last=.true. do k=1,imjm if (isea(k) .ne.0. .and. isea(k).ne. 1.) then write(6,*) 'trying to fix k= ', k last=.false. if (k .le. im) then !Bottom row if (k .eq. 1) then !Lower-left corner isea(k)=isea(k+im) elseif (k .eq. im) then !Lower-right corner isea(k)=isea(k+im1) else isea(k)=max(isea(k+im1),isea(k+im)) endif elseif (k .gt. imjm-im) then !Top row if (k .eq. imjm-im1) then !Upper-left corner isea(k)=isea(k-im1) elseif (k .eq. imjm) then !Upper-right corner isea(k)=isea(k-im) else isea(k)=max(isea(k-im1),isea(k-im)) endif else !Interior row if (mod(k,imt) .eq. 1) then !Left column, odd row isea(k)=max(isea(k-im1),isea(k+1),isea(k+im)) elseif (mod(k,imt) .eq. im) then !Right column, odd row isea(k)=max(isea(k-im),isea(k-1),isea(k+im1)) else isea(k)=max(isea(k+im1),isea(k+im) . ,isea(k-im1),isea(k-im)) endif endif endif enddo enddo Cmptest c c write(6,*) 'GOT TO THIS POINT' C open(1,file=soil_out,status='unknown',form='unformatted') C write(1) isea C close(1) call conh12t(isea,imjm,1,dum2d,imt,jmt) do J=JMT,1,-2 C do J=1,JMT,2 write(6,273)(dum2d(I,J),I=1,IMT,IMT/29) enddo 273 format(31f2.0) c return end C call corners(im,jm,imjm,tph0d,tlm0d,dlmd,dphd,apmn,almn,apmx,almx) subroutine CORNERS(im,jm,IMJM,tph0d,tlm0d,dlmd,dphd,minlat,wlon, + maxlat,elon) C C C * ROUTINE TO FIND EARTH LATITUDE/LONGITUDE FOR THE CORNER * C * POINTS OF AN ETA MODEL GRID * C C----------------------------------------------------------------------- D I M E N S I O N & KHL0 (JM),KHH0 (JM), GLAT(IMJM),GLON(IMJM) D A T A & PI/3.141592654/ REAL DLMD,DPHD,WBD,SBD,TPH0D,TLM0D,minlat,wlon,maxlat,elon C******************************* WBD=-(float(IM)-1.)*DLMD SBD=(-(float(JM)-1.)/2.)*DPHD DTR=PI/180. TPH0=TPH0D*DTR WB=WBD*DTR SB=SBD*DTR DLM=DLMD*DTR DPH=DPHD*DTR TDLM=DLM+DLM TDPH=DPH+DPH C STPH0=SIN(TPH0) CTPH0=COS(TPH0) C DO 100 J=1,JM KHL0(J)=IM*(J-1)-(J-1)/2+1 KHH0(J)=IM*J-J/2 C WRITE(6,9999) J, KHL0(J), KHH0(J) C9999 FORMAT(2X,3(I10,1X)) 100 CONTINUE C--------------GEOGRAPHIC LAT AND LONG OF TLL GRID POINTS--------------- TPH=SB-DPH maxlat=-999. minlat=99. DO 200 J=1,JM KHL=KHL0(J) KHH=KHH0(J) C TLM=WB-TDLM+MOD(J+1,2)*DLM TPH=TPH+DPH STPH=SIN(TPH) CTPH=COS(TPH) DO 200 K=KHL,KHH TLM=TLM+TDLM SPH=CTPH0*STPH+STPH0*CTPH*COS(TLM) GLAT(K)=ASIN(SPH) CLM=CTPH*COS(TLM)/(COS(GLAT(K))*CTPH0)-TAN(GLAT(K))*TAN(TPH0) IF(CLM.GT.1.) CLM=1. FACT=1. IF(TLM.GT.0.) FACT=-1. GLON(K)=(-TLM0D*DTR+FACT*ACOS(CLM))/DTR Cmp at this point GLON is in DEGREES WEST if (GLON(K) .lt. 0) GLON(K)=GLON(K)+360. if (GLON(K) .gt. 360.) GLON(K)=GLON(K)-360. if (GLON(K) .lt. 180) GLON(K)=-GLON(K) ! make WH negative if (GLON(K) .gt. 180) GLON(K)=360.-GLON(K) ! make EH GLAT(K)=GLAT(K)/DTR if (glat(k) .gt. maxlat) maxlat=glat(k) if (glat(k) .lt. minlat) minlat=glat(k) 200 CONTINUE if (TPH0D .ge. 0) then wlon=glon(imjm-im+1) elon=glon(imjm) else wlon=glon(1) elon=glon(im) endif C write(6,*) 'raw lon values (w,e) ', wlon, elon if (tlm0d .lt. 0 .and. wlon .gt. 0) wlon=wlon-360. if (tlm0d .gt. 0 .and. elon .lt. 0) elon=elon+360. RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SWAPVAL(A) C C REVERSE ORDER OF BYTES IN INTEGER*2 WORD, or REAL*2 C INTEGER*2 A C CHARACTER*1 JTEMP(2) CHARACTER*1 KTEMP C EQUIVALENCE (JTEMP(1),ITEMP) C ITEMP = A KTEMP = JTEMP(2) JTEMP(2) = JTEMP(1) JTEMP(1) = KTEMP A = ITEMP RETURN END C add subroutines subroutine conh12t(h1,imjm,lm,h2,imt,jmt) c c *** Routine for reordering thibue-height-point-1-dimensional c matrices for 2-dimensional indexing (imt,jmt). c implicit none c integer*4 imt,jmt,imjm,lm,i,j,k,l c real*4 h1(imjm,lm),h2(imt,jmt,lm) c_______________________________________________________________________________ c do l=1,lm k=0 do j=1,jmt do i=1,imt,2 k=k+1 h2(i,j,l)=h1(k,l) enddo c if (j .lt. jmt) then do i=2,imt,2 k=k+1 h2(i,j,l)=h1(k,l) enddo endif enddo enddo c return end c c=============================================================================== c=============================================================================== c subroutine conh12(h1,imjm,lm,h2,im,jm) c c *** Routine for reordering thibue-height-point-1-dimensional c matrices (im,jm) for 2-dimensional indexing. c implicit none c integer*4 im,jm,imjm,lm,i,j,k,l,imm1 c real*4 h1(imjm,lm),h2(im,jm,lm) c_______________________________________________________________________________ c imm1=im-1 do l=1,lm k=0 do j=1,jm do i=1,imm1+mod(j,2) k=k+1 h2(i,j,l)=h1(k,l) enddo enddo do j=2,jm-1,2 h2(im,j,l)=h2(imm1,j,l) enddo enddo c return end c c===============================================================================