! Copyright 2014 College of William and Mary ! ! Licensed under the Apache License, Version 2.0 (the "License"); ! you may not use this file except in compliance with the License. ! You may obtain a copy of the License at ! ! http://www.apache.org/licenses/LICENSE-2.0 ! ! Unless required by applicable law or agreed to in writing, software ! distributed under the License is distributed on an "AS IS" BASIS, ! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. ! See the License for the specific language governing permissions and ! limitations under the License. !Modified to output time stamp as 'yyyymmddhhmmss', and .csv format for Arc (2D !variables only) !**************************************************************************************** ! * ! read_output7b with (x,y,z) read in from station.bp or station.sta (sta format; x,y, with z), ! (z>=0 is distance from F.S.; if below bottom, the bottom value is returned) ! for 3D variables (surface values for 2D variables). ! Outputs: fort.1[89]; ; fort.20 - local dapth for each pt. ! For ics=2 (lat/lon), use nearest node for output ! * ! ifort -Bstatic -assume byterecl -O3 -o read_output7b_group_zfs2 read_output7b_group_zfs2.f90 ! PGI compiler: ! pgf90 -O2 -mcmodel=medium -Bstatic -o read_output7b_group_zfs2 read_output7b_group_zfs2.f90 ! On Ranch: ! f90 -Bdynamic -o read_output7b_group_zfs2 read_output7b_group_zfs2.f90 !**************************************************************************************** ! program read_out parameter(nbyte=4) ! parameter(mnp=80000) ! parameter(mne=160000) ! parameter(mnv=100) character*30 file63 character*12 it_char ! character*14 timestamp character*48 start_time,version,variable_nm,variable_dim character*48 data_format character*2 jday_char,jhour_char,jmin_char,jsec_char character(len=5) :: i_char(6000) integer(kind=8) :: itimestamp integer,allocatable :: elnode(:,:) ! dimension sigma(mnv),cs(mnv),ztot(mnv),x(mnp),y(mnp),dp(mnp),kbp00(mnp),kfp(mnp) ! dimension elnode(mne,4),out(mnp,3,mnv,2),out2(mnp,mnv,2),icum(mnp,mnv,2),eta2(mnp),node3(mnp,3),arco(mnp,3) ! dimension ztmp(mnv),x00(mnp),y00(mnp),iep(mnp),swild(3),out3(mnp,2),z00(mnp) dimension swild(3) allocatable :: sigma(:),cs(:),ztot(:),x(:),y(:),dp(:),kbp00(:),kfp(:) allocatable :: out(:,:,:,:),out2(:,:,:),icum(:,:,:),eta2(:),node3(:,:),arco(:,:) allocatable :: ztmp(:),x00(:),y00(:),iep(:),out3(:,:),z00(:),rl2min(:),dep(:) print*, 'Input extraction pts format (1: .bp; 2:.sta):' read(*,*)ibp if(ibp/=1.and.ibp/=2) stop 'Unknown format' print*, 'Input ics (1-Cartesian; 2-lat/lon):' read(*,*)ics print*, 'Input file to read from (without *_):' read(*,'(a30)')file63 print*, 'Input start and end file # to read:' read(*,*) iday1,iday2 print*, 'Input start CORIE day (0 if none):' read(*,*) corieday if(ibp==1) then !.bp format open(10,file='station.bp',status='old') read(10,*) read(10,*) nxy else !.sta format open(10,file='station.sta',status='old') read(10,*) nxy endif !ibp allocate(x00(nxy),y00(nxy),z00(nxy),rl2min(nxy),dep(nxy),stat=istat) if(istat/=0) stop 'Falied to allocate (1)' do i=1,nxy if(ibp==1) then !.bp format read(10,*)j,x00(i),y00(i),z00(i) else !.sta format read(10,*) read(10,*)x00(i),y00(i),z00(i) endif !ibp if(z00(i)<0) then write(*,*)'Invalid z value:',i; stop endif enddo !i close(10) !... Header !... open(63,file='1_'//file63,status='old',access='direct',recl=nbyte) irec=0 do m=1,48/nbyte read(63,rec=irec+m) data_format(nbyte*(m-1)+1:nbyte*m) enddo if(data_format.ne.'DataFormat v5.0') then print*, 'This code reads only v5.0: ',data_format stop endif irec=irec+48/nbyte do m=1,48/nbyte read(63,rec=irec+m) version(nbyte*(m-1)+1:nbyte*m) enddo irec=irec+48/nbyte do m=1,48/nbyte read(63,rec=irec+m) start_time(nbyte*(m-1)+1:nbyte*m) enddo irec=irec+48/nbyte do m=1,48/nbyte read(63,rec=irec+m) variable_nm(nbyte*(m-1)+1:nbyte*m) enddo irec=irec+48/nbyte do m=1,48/nbyte read(63,rec=irec+m) variable_dim(nbyte*(m-1)+1:nbyte*m) enddo irec=irec+48/nbyte write(*,'(a48)')data_format write(*,'(a48)')version write(*,'(a48)')start_time write(*,'(a48)')variable_nm write(*,'(a48)')variable_dim read(63,rec=irec+1) nrec read(63,rec=irec+2) dtout read(63,rec=irec+3) nspool read(63,rec=irec+4) ivs read(63,rec=irec+5) i23d irec=irec+5 print*, 'i23d=',i23d,' nrec= ',nrec ! Vertical grid read(63,rec=irec+1) nvrt read(63,rec=irec+2) kz read(63,rec=irec+3) h0 read(63,rec=irec+4) h_s read(63,rec=irec+5) h_c read(63,rec=irec+6) theta_b read(63,rec=irec+7) theta_f irec=irec+7 allocate(ztot(nvrt),sigma(nvrt),cs(nvrt),ztmp(nvrt),stat=istat) if(istat/=0) stop 'Falied to allocate (2)' do k=1,kz-1 read(63,rec=irec+k) ztot(k) enddo do k=kz,nvrt kin=k-kz+1 read(63,rec=irec+k) sigma(kin) cs(kin)=(1-theta_b)*sinh(theta_f*sigma(kin))/sinh(theta_f)+ & &theta_b*(tanh(theta_f*(sigma(kin)+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5) enddo irec=irec+nvrt ! Horizontal grid read(63,rec=irec+1) np read(63,rec=irec+2) ne irec=irec+2 allocate(x(np),y(np),dp(np),kbp00(np),kfp(np),elnode(4,ne),out(np,3,nvrt,2), & &out2(np,nvrt,2),icum(np,nvrt,2),eta2(np),node3(np,3),arco(np,3), & &iep(np),out3(np,2),stat=istat) if(istat/=0) stop 'Falied to allocate (3)' do m=1,np read(63,rec=irec+1)x(m) read(63,rec=irec+2)y(m) read(63,rec=irec+3)dp(m) read(63,rec=irec+4)kbp00(m) irec=irec+4 enddo !m=1,np do m=1,ne read(63,rec=irec+1)i34 irec=irec+1 do mm=1,i34 read(63,rec=irec+1)elnode(mm,m) irec=irec+1 enddo !mm enddo !m irec0=irec print*, 'last element',(elnode(j,ne),j=1,3) !... Find parent element for (x00,y00) iep=0 if(ics==1) then !Cartesian do i=1,ne do l=1,nxy aa=0 ar=0 !area do j=1,3 j1=j+1 j2=j+2 if(j1>3) j1=j1-3 if(j2>3) j2=j2-3 n0=elnode(j,i) n1=elnode(j1,i) n2=elnode(j2,i) swild(j)=signa(x(n1),x(n2),x00(l),y(n1),y(n2),y00(l)) !temporary storage aa=aa+abs(swild(j)) if(j==1) then ar=signa(x(n1),x(n2),x(n0),y(n1),y(n2),y(n0)) if(ar<=0) then print*, 'Negative area:',ar,i,n1,n2,n0,x(n1),y(n1),x(n2),y(n2),x(n0),y(n0) stop endif endif enddo !j ae=abs(aa-ar)/ar if(ae<=1.e-5) then iep(l)=i node3(l,1:3)=elnode(1:3,i) arco(l,1:3)=swild(1:3)/ar arco(l,1)=max(0.,min(1.,arco(l,1))) arco(l,2)=max(0.,min(1.,arco(l,2))) if(arco(l,1)+arco(l,2)>1) then arco(l,3)=0 arco(l,2)=1-arco(l,1) else arco(l,3)=1-arco(l,1)-arco(l,2) endif cycle endif enddo !l; build pts ifl=0 !flag do l=1,nxy if(iep(l)==0) then ifl=1 exit endif enddo !l if(ifl==0) exit enddo !i=1,ne else !lat/lon rl2min=1.e25 !min distance^2 do i=1,np do l=1,nxy rl2=(x(i)-x00(l))**2+(y(i)-y00(l))**2 if(rl26000) stop 'Increase dim' do i=1,6000 write(i_char(i),'(i5.5)')i enddo !i write(18,'(a1,6000(a5,a5))')'t',(',elev',i_char(i),i=1,nxy) if(ivs==2) write(19,'(a1,6000(a5,a5))')'t',(',elev',i_char(i),i=1,nxy) !... Time iteration !... it_tot=(iday1-1)*nrec do iday=iday1,iday2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ write(it_char,'(i12)')iday it_char=adjustl(it_char) leng=len_trim(it_char) open(63,file=it_char(1:leng)//'_'//file63,status='old',access='direct',recl=nbyte) !' irec=irec0 do it1=1,nrec read(63,rec=irec+1) time read(63,rec=irec+2) it irec=irec+2 it_tot=it_tot+1 time=it_tot*dtout ! Timestamp if(time>86400) stop 'Exceed a day' ! jday=time/86400 jhour=time/3600 jmin=(time-jhour*3600)/60 jsec=time-jhour*3600-jmin*60 ! write(jday_char,'(i2.2)')jday+1 ! write(jhour_char,'(i2.2)')jhour ! write(jmin_char,'(i2.2)')jmin ! write(jsec_char,'(i2.2)')jsec !timestamp='20010101'//jhour_char//jmin_char//jsec_char itimestamp=20010101000000_8+jhour*10000+jmin*100+jsec ! print*, 'time=',time/86400 do i=1,nxy do j=1,3 nd=node3(i,j) read(63,rec=irec+nd) eta2(nd) enddo !j enddo !i ! do i=1,np ! read(63,rec=irec+i) eta2(i) ! enddo !i irec=irec+np out2=0 out3=0 if(i23d.eq.2) then do i=1,nxy dep(i)=0 do j=1,3 !nodes nd=node3(i,j) !Compute local depth dep(i)=dep(i)+arco(i,j)*dp(nd) do m=1,ivs read(63,rec=irec+(nd-1)*ivs+m) tmp out2(i,1,m)=out2(i,1,m)+arco(i,j)*tmp enddo !m enddo !j enddo !i irec=irec+np*ivs write(18,'(i14,6000(a1,e14.6))')itimestamp,(',',out2(i,1,1),i=1,nxy) if(ivs==2) write(19,'(i14,6000(a1,e14.6))')itimestamp,(',',out2(i,1,2),i=1,nxy) else !i23d=3 do i=1,nxy do j=1,3 !nodes nd=node3(i,j) do k=max0(1,kbp00(nd)),nvrt do m=1,ivs read(63,rec=irec+icum(nd,k,m)) out(i,j,k,m) enddo !m enddo !k enddo !j enddo !i irec=irec+icum(np,nvrt,ivs) ! Do interpolation do i=1,nxy etal=0; dep(i)=0; idry=0 do j=1,3 nd=node3(i,j) if(eta2(nd)+dp(nd)=ztot(k).and.-dep(i)=ztmp(k).and.ztmp(nvrt)-z00(i)<=ztmp(k+1)) then k0=k rat=(ztmp(nvrt)-z00(i)-ztmp(k))/(ztmp(k+1)-ztmp(k)) exit endif enddo !k endif !ztmp if(k0==0) then ! write(12,*)'Warning: failed to find a vertical level:',it,i else do m=1,ivs out3(i,m)=out2(i,k0,m)*(1-rat)+out2(i,k0+1,m)*rat enddo !m endif endif !dry/wet enddo !i=1,nxy write(18,'(e16.8,6000(1x,f12.3))')time/86400+corieday,(out3(i,1),i=1,nxy) if(ivs==2) write(19,'(e16.8,6000(1x,f12.3))')time/86400+corieday,(out3(i,2),i=1,nxy) endif !i23d enddo !it1=1,nrec !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ enddo !iday ! Output local depths info do i=1,nxy write(20,*)i,dep(i) enddo !i stop end ! function signa(x1,x2,x3,y1,y2,y3) !!... Compute signed area formed by pts 1,2,3 ! ! signa=((x1-x3)*(y2-y3)-(x2-x3)*(y1-y3))/2 ! ! return ! end