program lagrtracer c The program calculates the moving of tracers with water using current velocity c and diffusing randomly the tracers depending on diffusion coefficient c !uses fortran random function specific for ifort, g77 or other, check this! c !to be included the wind stress effects c assumes: when outside model area- put undef number -9999. c assumes: when reache the coast - stay there c c As input: c This program reads as input parameters: oilspill_parameters.inp c number of tracers to realease continuously in time c initial point of release of tracer, also the spot size c tracers move with the water parcels c model simulation - field of U and V component of currents c As output: c files named fort.40x for every moment+1 of location of tracers - lon;lat c For example 10 moments of time: 11 files fort.401 fort.402 ... fort.411 parameter (pi=3.1415926536,re=6378.1e3,undef=-9999.) real,allocatable :: u(:,:,:,:),v(:,:,:,:) real,allocatable :: alon(:),alat(:),time(:),mask(:,:) real,allocatable :: tlon(:,:),tlat(:,:) integer rc,ncid,ndims,imt,jmt,km,nt integer lat_id,lon_id,time_id,u_id,v_id integer dims(4),start4d(4),edges4d(4) integer iseed/2011/ character*10 name_u,name_v character*500 ifile include "netcdf.inc" name_u='uu' !eastward_sea_water_velocity name_v='vv' !northward_sea_water_velocity include "oilspill_parameters.inp" c lt - number of lagrangian tracers to let move c ah - horizontal diffusion coefficient ~10-100 m2/s c the initial coordinates of the point of the oil spill c the initial radius of the oil spill in meters, then convert in degrees c ifile - nc file to read the surface current velocity from [m/s] c !wind stress not yet! to be included!!! c oillon=9.7 c oillat=40.25 c radstart=1000. c write (*,*) 'input file to read currents =?' c read (*,*) ifile c write (*,*) 'ifile=',ifile radstartx=radstart/(pi/180.*re*cos(oillat*pi/180.)) radstarty=radstart/(pi/180.*re) raddeg=sqrt(radstartx*radstartx+radstarty*radstarty) c read the model grid, velocity field and make mask err=nf_open(ifile,NCNOWRIT,ncid) if (err .ne. NF_NOERR) write (*,*) 'error in opening ',ifile err = nf_inq_ndims(ncid,ndims) if (err .NE. NF_NOERR) write (*,*) 'error in ndims ' err = nf_inq_dimlen(ncid,1,imt) if (err .ne. NF_NOERR) write (*,*) 'error in imt' err = nf_inq_dimlen(ncid,2,jmt) if (err .ne. NF_NOERR) write (*,*) 'error in jmt' err = nf_inq_dimlen(ncid,3,km) if (err .ne. NF_NOERR) write (*,*) 'error in km' err = nf_inq_dimlen(ncid,4,nt) if (err .ne. NF_NOERR) write (*,*) 'error in nt ' write (*,*) 'imt,jmt,km,nt',imt,jmt,km,nt start4d(1)=1 edges4d(1) = imt start4d(2)=1 edges4d(2) = jmt start4d(3)=1 edges4d(3) = km start4d(4)=1 edges4d(4) = nt c read longitude allocate (alon(imt),stat=rc) if (rc /= 0) write (*,*) ' Error allocating alon' err = nf_inq_varid(ncid,'lonc',lon_id) if (err .NE. NF_NOERR) write (*,*) 'error in lon_id' err = nf_get_vara_real(ncid,lon_id,1,imt,alon) if (err .ne. NF_NOERR) write (*,*) 'error in get alon data' c read latitude allocate (alat(jmt),stat=rc) if (rc /= 0) write (*,*) ' Error allocating alat' err = nf_inq_varid(ncid,'latc',lat_id) if (err .NE. NF_NOERR) write (*,*) 'error in lat_id' err = nf_get_vara_real(ncid,lat_id,1,jmt,alat) if (err .ne. NF_NOERR) write (*,*) 'error in get alat data' c read time allocate (time(nt),stat=rc) if (rc /= 0) write (*,*) ' Error allocating time' err = nf_inq_varid(ncid,'time',time_id) if (err .NE. NF_NOERR) write (*,*) 'error in time_id' err = nf_get_vara_real(ncid,time_id,1,nt,time) if (err .ne. NF_NOERR) write (*,*) 'error in get time data' c read u values err = nf_inq_varid(ncid,name_u,u_id) if (err .NE. NF_NOERR) write (*,*) 'error in u_id' allocate (u(imt,jmt,km,nt),stat=rc) if (rc /= 0) write (*,*) ' Error allocating u' err = nf_get_vara_real(ncid,u_id,start4d,edges4d,u) if (err .ne. NF_NOERR) write (*,*) 'error in get u data' c read v values err = nf_inq_varid(ncid,name_v,v_id) if (err .NE. NF_NOERR) write (*,*) 'error in v_id' allocate (v(imt,jmt,km,nt),stat=rc) if (rc /= 0) write (*,*) ' Error allocating v' err = nf_get_vara_real(ncid,v_id,start4d,edges4d,v) if (err .ne. NF_NOERR) write (*,*) 'error in get v data' c close the input file err=nf_close(ncid) if (err .ne. NF_NOERR) write (*,*) 'error in closing ',ifile allocate (mask(imt,jmt),stat=rc) if (rc /= 0) write (*,*) ' Error allocating mask' allocate (tlon(nt+1,lt*nt),stat=rc) if (rc /= 0) write (*,*) ' Error allocating tlon' allocate (tlat(nt+1,lt*nt),stat=rc) if (rc /= 0) write (*,*) ' Error allocating tlat' c create the mask do i=1,imt do j=1,jmt if (u(i,j,km,1).ne.undef.and.v(i,j,km,1).ne.undef) then mask(i,j)=1. else mask(i,j)=0. endif enddo enddo startlon=alon(1) startlat=alat(1) endlon=alon(imt) endlat=alat(jmt) deltat=time(2)-time(1) write (*,*) 'deltal=',deltat c calculate the factor for diffusion diff=sqrt(6.*ah*deltat) write (*,*) 'Area between:',startlon, endlon,' E and + ',startlat, endlat,' N' write (*,*) 'Oil spill location:',oillon,oillat c calculate the initial positions of tracers using random number call random_seed(size=iseed) do l=1,lt c for f77 function rand c tlon(1,l)=oillon+(2.*rand(l)-1.)*raddeg* c + asin(radstarty/radstartx) c tlat(1,l)=oillat+(2.*rand(l)-1.)*raddeg* c + acos(radstarty/radstartx) c for ifort function ran c tlon(1,l)=oillon+(2.*ran(iseed)-1.)*raddeg* c + asin(radstarty/radstartx) c tlat(1,l)=oillat+(2.*ran(iseed)-1.)*raddeg* c + acos(radstarty/radstartx) c for fortran 90 (pgf90, gfortran) function random_number call random_number(harvest=rannum) tlon(1,l)=oillon+(2.*rannum-1.)*raddeg* + asin(radstarty/radstartx) call random_number(harvest=rannum) tlat(1,l)=oillat+(2.*rannum-1.)*raddeg* + acos(radstarty/radstartx) if (tlon(1,l).lt.startlon.or.tlon(1,l).gt.endlon.or. + tlat(1,l).lt.startlat.or.tlat(1,l).gt.endlat) then tlon(1,l)=undef tlat(1,l)=undef endif enddo c calculate the next positions of tracers do 20 n=2,nt+1 do 10 l=1,(n-1)*lt if (tlon(n-1,l).eq.undef) then tlon(n,l)=undef tlat(n,l)=undef goto 10 endif do i=1,imt-1 do j=1,jmt-1 if (tlon(n-1,l).ge.alon(i).and.tlon(n-1,l).lt.alon(i+1)) then if (tlat(n-1,l).ge.alat(j).and.tlat(n-1,l).lt.alat(j+1)) then if (mask(i,j)+mask(i+1,j)+mask(i,j+1)+mask(i+1,j+1).ne.4) then tlon(n,l)=tlon(n-1,l) tlat(n,l)=tlat(n-1,l) goto 10 else c calculate the velocity in the point w1=(tlon(n-1,l)-alon(i))*cos(tlat(n-1,l)*pi/180.)* + (tlat(n-1,l)-alat(j)) w2=(tlon(n-1,l)-alon(i))*cos(tlat(n-1,l)*pi/180.)* + (alat(j+1)-tlat(n-1,l)) w3=(alon(i+1)-tlon(n-1,l))*cos(tlat(n-1,l)*pi/180.)* + (tlat(n-1,l)-alat(j)) w4=(alon(i+1)-tlon(n-1,l))*cos(tlat(n-1,l)*pi/180.)* + (alat(j+1)-tlat(n-1,l)) w=w1+w2+w3+w4 tu=1./w*(w4*u(i,j,km,n-1)+w2*u(i+1,j,km,n-1)+ + w1*u(i+1,j+1,km,n-1)+w3*u(i,j+1,km,n-1)) tv=1./w*(w4*v(i,j,km,n-1)+w2*v(i+1,j,km,n-1)+ + w1*v(i+1,j+1,km,n-1)+w3*v(i,j+1,km,n-1)) c advection tlon(n,l)=tlon(n-1,l)+tu*deltat/ + (pi/180.*re*cos(tlat(n-1,l)*pi/180.)) tlat(n,l)=tlat(n-1,l)+tv*deltat/(pi/180.*re) c if (l.eq.1) write (*,*) tlon(n,l),tlat(n,l) c diffusion call random_number(harvest=rannum) tlon(n,l)=tlon(n,l)+(2.*rannum-1.)*diff/ + (pi/180.*re*cos(tlat(n-1,l)*pi/180.)) call random_number(harvest=rannum) tlat(n,l)=tlat(n,l)+(2.*rannum-1.)*diff/(pi/180.*re) c wind stress c ??? goto 15 endif endif endif enddo enddo 15 if (tlon(n,l).lt.startlon.or.tlon(n,l).gt.endlon.or. + tlat(n,l).lt.startlat.or.tlat(n,l).gt.endlat) then tlon(n,l)=undef tlat(n,l)=undef endif 10 continue l0=1 do l=(n-1)*lt+1,n*lt tlon(n,l)=tlon(1,l0) tlat(n,l)=tlat(1,l0) l0=l0+1 enddo 20 continue c write results c one file with the position of each tracer, commented for now c do l=1,lt c do n=1,nt+1 c if (tlon(n,l).ne.undef) write (100+l,*) tlon(n,l),tlat(n,l) c enddo c enddo c one file for every moment of time, plus the next one, e.g. nt+1 files do n=1,nt+1 do l=1,n*lt if (tlon(n,l).ne.undef) write (400+n,*) tlon(n,l),tlat(n,l) enddo enddo deallocate (alon,stat=rc) if (rc /= 0) write (*,*) ' Error deallocating alon' deallocate (alat,stat=rc) if (rc /= 0) write (*,*) ' Error deallocating alat' deallocate (tlon,stat=rc) if (rc /= 0) write (*,*) ' Error deallocating tlon' deallocate (tlat,stat=rc) if (rc /= 0) write (*,*) ' Error deallocating tlat' deallocate (u,stat=rc) if (rc /= 0) write (*,*) ' Error deallocating u' deallocate (v,stat=rc) if (rc /= 0) write (*,*) ' Error deallocating v' deallocate (mask,stat=rc) if (rc /= 0) write (*,*) ' Error deallocating mask' stop end