program extract_HC ccc LANA, 2004 (based on old heval and uveval) ccc reads OTIS standard binary complex model file ccc (elevations OR transports), reads a list of locations ccc and outputs ASCII file with the complex amplitudes/amp,phases of ccc elevations/transports/currents interpolated at the locations c implicit none include 'constit.h' complex, allocatable:: z(:,:,:),z1(:),dtmp(:,:) complex, allocatable:: zl(:,:,:),zl1(:) complex d1 real, allocatable:: lat(:),lon(:),depth(:,:),x(:),y(:),lon0(:) real latp,lonp,xt,yt real th_lim(2),ph_lim(2),dum,lth_lim(2),lph_lim(2) integer, allocatable:: mask(:,:),cind(:),lcind(:), * mz(:,:),mzl(:,:) c character*4 c_id(ncmx),c_id_mod(ncmx),lc_id(ncmx) character*80 modname,lltname,outname,gname,ctmp,lname character*2000 fmt character*80 rmCom character*1 zuv,c1 character*80 xy_ll_sub logical APRI,geo,interp_micon integer ncon,nc,n,m,ndat,k,ierr,ierr1,ic,n0,m0 integer ncl,nl,ml,nmod,imod,ibl c nmod=1 ibl=0 lname='DATA/load_file' call rd_inp(modname,lltname,zuv,c_id,ncon,APRI,geo, * outname,interp_micon) if(trim(modname).eq.'model.list')then nmod=0 open(unit=12,file='model.list',status='old',err=11) write(*,*)'Models to extract HC from:' 13 read(12,'(a80)',end=12)ctmp write(*,*)trim(ctmp) nmod=nmod+1 go to 13 12 rewind(12) read(12,'(a80)')modname endif write(*,*) write(*,*)'Lat/Lon file:',trim(lltname) if(ncon.gt.0)write(*,*)'Constituents to include: ',c_id(1:ncon) if(geo.and.zuv.eq.'z')then write(*,*)'Extract GEOCENTRIC tide HC' else write(*,*)'Extract OCEAN tide HC' endif c open(unit=11,file=outname,status='unknown') do imod=1,nmod write(*,*) if(imod.gt.1)read(12,'(a80)')modname call rd_mod_header(modname,zuv,n,m,th_lim,ph_lim,nc,c_id_mod, * xy_ll_sub) write(*,*)'Model: ',trim(modname(12:80)) write(11,*)'Model: ',trim(modname(12:80)) if(trim(xy_ll_sub).eq.'')then write(*,*)'Lat limits: ',th_lim write(*,*)'Lon limits: ',ph_lim else if(trim(xy_ll_sub).ne.'xy_ll_N'.and. * trim(xy_ll_sub).ne.'xy_ll_S')then write(*,*)'No converting function ', trim(xy_ll_sub), * 'in the OTPS' stop endif if(trim(xy_ll_sub).eq.'xy_ll_N')then call xy_ll_N(ph_lim(1),th_lim(1),xt,yt) if(xt.gt.180)xt=xt-360 write(*,*)'Lat,Lon lower left corner:',yt,xt call xy_ll_N(ph_lim(1),th_lim(2),xt,yt) if(xt.gt.180)xt=xt-360 write(*,*)'Lat,Lon upper left corner:',yt,xt call xy_ll_N(ph_lim(2),th_lim(1),xt,yt) if(xt.gt.180)xt=xt-360 write(*,*)'Lat,Lon lower right corner:',yt,xt call xy_ll_N(ph_lim(2),th_lim(2),xt,yt) if(xt.gt.180)xt=xt-360 write(*,*)'Lat,Lon upper right corner:',yt,xt else call xy_ll_S(ph_lim(1),th_lim(1),xt,yt) if(xt.gt.180)xt=xt-360 write(*,*)'Lat,Lon lower left corner:',yt,xt call xy_ll_S(ph_lim(1),th_lim(2),xt,yt) if(xt.gt.180)xt=xt-360 write(*,*)'Lat,Lon upper left corner:',yt,xt call xy_ll_S(ph_lim(2),th_lim(1),xt,yt) if(xt.gt.180)xt=xt-360 write(*,*)'Lat,Lon lower right corner:',yt,xt call xy_ll_S(ph_lim(2),th_lim(2),xt,yt) if(xt.gt.180)xt=xt-360 write(*,*)'Lat,Lon upper right corner:',yt,xt endif endif write(*,*)'Constituents: ',c_id_mod(1:nc) if(trim(xy_ll_sub).ne.'')then write(*,*)'Model is on uniform grid in km' write(*,*)'Function to convert x,y to lat,lon:', * trim(xy_ll_sub) endif c if(zuv.eq.'z')write(*,*)'Output elevations (m)' if(zuv.eq.'U')write(*,*)'Output WE transport (m^2/s)' if(zuv.eq.'u')write(*,*)'Output WE velocity (cm/s)' if(zuv.eq.'V')write(*,*)'Output SN transport (m^2/s)' if(zuv.eq.'v')write(*,*)'Output SN velocity (cm/s)' c if(zuv.eq.'z')write(11,*)'Elevations (m)' if(zuv.eq.'U')write(11,*)'WE transport (m^2/s)' if(zuv.eq.'u')write(11,*)'WE velocity (cm/s)' if(zuv.eq.'V')write(11,*)'SN transport (m^2/s)' if(zuv.eq.'v')write(11,*)'SN velocity (cm/s)' c if(ncon.eq.0)then ibl=1 ncon=nc c_id=c_id_mod write(*,*)'Constituents to include: ',c_id(1:ncon) endif c allocate(cind(ncon)) call def_con_ind(c_id,ncon,c_id_mod,nc,cind) c ndat=0 open(unit=1,file=lltname,status='old',err=1) 3 read(1,*,end=2)dum,dum ndat=ndat+1 go to 3 2 rewind(1) allocate(lat(ndat),lon(ndat),lon0(ndat)) if(trim(xy_ll_sub).ne.'')allocate(x(ndat),y(ndat)) do k=1,ndat read(1,*)lat(k),lon(k) if(trim(xy_ll_sub).eq.'xy_ll_N')then call ll_xy_N(lon(k),lat(k),x(k),y(k)) elseif(trim(xy_ll_sub).eq.'xy_ll_S')then call ll_xy_S(lon(k),lat(k),x(k),y(k)) endif lon0(k)=lon(k) if(trim(xy_ll_sub).eq.'')then ! check on lon convention if(lon(k).gt.ph_lim(2))lon(k)=lon(k)-360 if(lon(k).lt.ph_lim(1))lon(k)=lon(k)+360 endif enddo ! k close(1) c allocate(z(ncon,n,m),z1(ncon),mask(n,m)) write(*,'(a,$)')' Reading model...' call rd_mod_body(modname,zuv,cind,ncon,z,n,m,mask) write(*,*)'done' if(zuv.eq.'z'.and.geo)then call rd_mod_header1(lname,nl,ml,ncl,lth_lim,lph_lim,lc_id) allocate(lcind(ncon)) call def_con_ind(c_id,ncon,lc_id,ncl,lcind) allocate(zl(ncon,nl,ml),zl1(ncon),mzl(nl,ml)) write(*,*) write(*,'(a,$)')' Reading file to apply load corrections...' call rd_mod_body1(lname,'z',lcind,ncon,zl,nl,ml,mzl) write(*,*)'done' endif if(zuv.eq.'u'.or.zuv.eq.'v')then allocate(depth(n,m),dtmp(n,m),mz(n,m)) c currents case: need to read grid open(unit=1,file=modname,status='old') read(1,*)gname read(1,*)gname read(1,'(a80)')gname gname=rmCom(gname) close(1) open(unit=1,file=trim(gname),status='old', * form='unformatted',err=4) read(1)n0,m0 ! ignore the rest of record if(n0.ne.n.or.m0.ne.m)then write(*,*)'Wrong grid in ',trim(gname) write(*,*)'Grid size and model size are different!' stop endif read(1) ! pass iobc read(1) depth read(1) mz close(1) dtmp=depth deallocate(depth) endif c output format fmt='(f9.3,f9.3,' c1=',' do ic=1,ncon if(ic.eq.ncon)c1=')' if(APRI)then fmt=trim(fmt)//'f8.3,f8.1'//c1 else fmt=trim(fmt)//'f8.3,f8.3'//c1 endif enddo c if(APRI)then write(11,*)' Lat Lon ', * ((trim(c_id(ic)),'_amp ', * trim(c_id(ic)),'_ph '),ic=1,ncon) else write(11,*)' Lat Lon ', * ((trim(c_id(ic)),'_Re ', * trim(c_id(ic)),'_Im '),ic=1,ncon) endif c1=zuv ! since interp change zuv (U->u, V->v) latp=0. lonp=0. do k=1,ndat if(lat(k).ne.latp.or.lon(k).ne.lonp)then if(trim(xy_ll_sub).eq.'')then call interp(z,ncon,n,m,mask,th_lim,ph_lim, * lat(k),lon(k),z1,ierr,c1) else call interp(z,ncon,n,m,mask,th_lim,ph_lim, * y(k),x(k),z1,ierr,c1) endif if(ierr.eq.0)then if(zuv.eq.'u'.or.zuv.eq.'v')then if(trim(xy_ll_sub).eq.'')then call interp(dtmp,1,n,m,mz,th_lim,ph_lim, * lat(k),lon(k),d1,ierr1,'z') else call interp(dtmp,1,n,m,mz,th_lim,ph_lim, * y(k),x(k),d1,ierr1,'z') endif z1=z1/real(d1)*100. ! currents cm/s elseif(zuv.eq.'z'.and.geo)then call interp(zl,ncon,nl,ml,mzl,lth_lim,lph_lim, * lat(k),lon(k),zl1,ierr1,'z') z1=z1+zl1 ! apply load correction to get geocentric tide endif endif if(ierr.eq.0)then if(APRI)then write(11,fmt)lat(k),lon0(k),((abs(z1(ic)),atan2(-imag(z1(ic)), * real(z1(ic)))*180/3.141593),ic=1,ncon) else write(11,fmt)lat(k),lon0(k), * ((real(z1(ic)),imag(z1(ic))),ic=1,ncon) endif else write(11,'(1x,f8.3,f8.3,a70)')lat(k),lon(k), * '************* Site is out of model grid OR land ***************' endif latp=lat(k) lonp=lon(k) endif enddo deallocate(z,z1,mask,cind,lat,lon,lon0) if(zuv.eq.'u'.or.zuv.eq.'v')deallocate(dtmp,mz) if(trim(xy_ll_sub).ne.'')deallocate(x,y) if(ibl.eq.1.and.geo)then ncon=0 deallocate(zl,zl1,mzl,lcind) endif enddo ! imod close(11) write(*,*)'Results are in ',trim(outname) if(geo.and.ibl.eq.0)deallocate(zl,zl1,mzl,lcind) stop 1 write(*,*)'Lat Lon file ',trim(lltname),' not found' write(*,*)'Check file setup.inp, line 2.' stop 4 write(*,*)'Grid file ',trim(gname),' not found' write(*,*)'Check file ',trim(modname) stop 11 write(*,*)'File ''model.list'' was NOT found...' write(*,*)'TO CREATE please do:' write(*,*)'ls -1 DATA/Model_*>model.list' stop end