module mod_floats use mod_xc ! HYCOM communication interface use mod_pipe ! HYCOM debugging interface !AS use mod_year_info, only : rt use mod_hycom_nersc, only : rungen #ifdef CALANUS ! For calanus model use mod_CAL06_calnusparticles use m_CAL06_initializecfmod use m_CAL06_read_restart use m_CAL06_calanus use mod_CAL06_read_norwecomcarbon #endif !AS end c implicit none c --- HYCOM synthetic floats, drifters and moorings c --- See subroutine floats (below) for more information c integer, parameter, public :: & nfldim=100001 !maximum number of synthetic floats real, allocatable, dimension(:,:,:), & save, public :: & wveli, ! interface vertical velocity & uold2, & vold2, & wold2u, & wold2d real, allocatable, dimension(:,:), & save, public :: & dlondx, & dlondy, & dlatdx, & dlatdy real, allocatable, dimension(:,:,:), & save, private :: & wvelup, & wveldn real, allocatable, dimension(:,:), & save, private :: & dpdxup, & dpdyup, & dpdxdn, & dpdydn !Varibles for the correction for curvliear grid real :: dlon, dlat, dlox, dlax, dloy, dlay, q_ll real, save :: flt(nfldim,13), & deltfl,fldepm,tbvar,tdecri,dtturb,uturb0, & mindepth(nfldim), maxdepth(nfldim) integer, save :: kfloat(nfldim),iflnum(nfldim),ifltyp(nfldim), & nflsam,nfldta,fltflg,nfladv,intpfl,iturbv,ismpfl, & nflout,nfl,nflt,nstepfl logical, save :: synflt,turbvel,samplfl,wvelfl,hadvfl,nonlatlon character*48 , save :: flnmflti,flnmflto,flnmfltio contains subroutine floats_init(m,n,time0) implicit none c include 'common_blocks.h' c integer m,n real time0 c c --- read initial float data c --- initialize parameters and arrays required for floats c integer i,j,k,l, nfl, new_nflt,itmp real rtmp integer nbcday,nsmday,ityp c include 'stmt_fns.h' c c --- initialize flags c c --- wvelfl: calculate vertical velocity? wvelfl=.false. c c --- hadvfl: horizontally advect the float? hadvfl=.false. c c --- nonlatlon: does the model grid cross latitude/longitude lines anywhere c --- in the domain? nonlatlon=.false. c c --- initialize file names flnmflti = 'floats.input' !AS flnmflto = 'floats_out' flnmfltio= 'floats.input_out' c c ------------------------- c --- set timing parameters c ------------------------- c c --- must have an integer number of baroclinic time steps per day nbcday=nint(86400.0/baclin) if (real(nbcday).ne.86400.0/baclin) then if (mnproc.eq.1) then write(lp,'(/ a /)') & 'error - need integer no. of bacl. time steps/day' call flush(lp) endif !1st tile call xcstop('(floats_init)') stop '(floats_init)' endif c c --- parameter nfladv is entered in blkdat.input as the number of c --- baroclinic time steps between updates of float position. at model c --- time nstep when the float is to be advected, velocity fields saved c --- at times nstep, nstep-nfladv/2, and nstep-nfladv are used to perform c --- the runga-kutta time interpolation. c c --- nfladv must be set so the float is advected every 2-4 hours, but c --- also must be no smaller than 4. This ensures that the runga- c --- kutta time interpolation is performed with a delta-time of 1-2 hours. c c --- if all floats are stationary (synthetic moorings), then nfladv c --- must still be no smaller than 4, but the above time restrection c --- is not necessary. this allows the user to sample at very high c --- frequency. for example, if baclin is 360 seconds, then setting c --- nfladv to 10 will allow water properties to be sampled once per hour c c --- set variable nfldta to nfladv/2 so that it equals the number of time c --- steps separating the velocity fields input into the runga-kutta c --- interpolation c nfldta=nfladv/2 c c --- make sure that the float will be advected at an integer number of c --- temporal points per day nsmday=nbcday/nfladv if (int(real(nbcday)/real(2*nfldta)).ne.nsmday) then if (mnproc.eq.1) then write(lp,'(/ a /)') & 'error - need integer no. of adv. times per day' call flush(lp) endif !1st tile call xcstop('(floats_init)') stop '(floats_init)' endif c c --- test to make sure that the advection time interval is between 2 c --- and 4 hours. only stop the program later if there are floats to c --- be advected if (nfladv.ne.4 .and. (nsmday.gt.12 .or. nsmday.lt.6)) then nfladv=-nfladv endif c c --- initialize number of time steps between float output if (nflsam.gt.0) then nflsam=nflsam*nbcday elseif (nflsam.eq.0) then nflsam=abs(nfladv) endif c c --- calculate the time interval for the runga-kutta interpolation c --- as 1/2 of the advection time interval deltfl=nfldta*baclin c c --- set minimum float depth (m) fldepm=0.0 c c --- ------------------------------------------------------------- c --- initialize synthetic drifters, floats, and moored instruments c --- ------------------------------------------------------------- c open(unit=uoff+99,file=trim(flnminp)//flnmflti,status='old') !on all nodes c c --- first line in float input file is the number of floats read(uoff+99,*) nflt if (nflt.gt.nfldim) then if (mnproc.eq.1) then write(lp,'(/ a /)') & 'error - increase nfldim in mod_floats.F' call flush(lp) endif !1st tile call xcstop('(floats_init)') stop '(floats_init)' endif if (nflt.eq.0) then synflt=.false. endif c c --- second line in float input file is the initial float time step c --- this number is set to zero for the first model run segment read(uoff+99,*) nstepfl c do nfl=1,nflt c c --- read input file containing the following information for each float: c --- c --- column 1 - initial sequential float number c --- column 2 - float type c --- 1 = 3-d lagrangian (vertically advected by diagnosed w) c --- 2 = isopycnic c --- 3 = isobaric (surface drifter when released in sfc. layer) c --- 4 = stationary (synthetic moored instrument) c --- column 3 - deployment time (days from model start, 0.0 = immediate) c --- column 4 - termination time (days from model start, 0.0 = forever) c --- column 5 - initial longitude (must be between minimum and maximum c --- longitudes defined in regional.grid.b) c --- column 6 - initial latitude (must be between minimum and maximum c --- latitudes defined in regional.grid.b) c --- column 7 - initial depth (or reference sigma for isopycnic floats) c read(uoff+99,*) iflnum(nfl),ifltyp(nfl),flt(nfl,8),flt(nfl,9), & flt(nfl,1),flt(nfl,2),flt(nfl,3) mindepth(nfl)=flt(nfl,3) maxdepth(nfl)=flt(nfl,3)+100.0; if (ifltyp(nfl).eq.1 .or. ifltyp(nfl).eq.4) then wvelfl=.true. endif if (ifltyp(nfl).ne.4) then hadvfl=.true. c c --- since this float is to be advected, stop if the advection time interval c --- is not between 2 and 4 hr if (nfladv.lt.0) then if (mnproc.eq.1) then write(lp,'(/ a /)') & 'error - set nfladv to advect floats every 2-4 hr' call flush(lp) endif !1st tile call xcstop('(floats_init)') stop '(floats_init)' endif c endif c flt(nfl,4)=0.0 flt(nfl,5)=0.0 flt(nfl,6)=0.0 flt(nfl,7)=0.0 if(ifltyp(nfl).eq.2) then flt(nfl,7)=flt(nfl,3) flt(nfl,3)=0.0 endif flt(nfl,10)=0.0 flt(nfl,11)=0.0 flt(nfl,12)=0.0 flt(nfl,13)=0.0 c kfloat(nfl)=-1 c enddo c nfladv=iabs(nfladv) c close(unit=uoff+99) !file='float.input' #if defined (CALANUS) if (mpop > nfldim) then if (mnproc.eq.1) then print*, 'mpop is larger than nfldim' call xcstop('(floats_init)') stop '(floats_init)' end if end if if (mnproc.eq.1) print*, 'initialize calanus, fltflg = ', fltflg if (mnproc.eq.1) then do nfl=1,nflt cfs(nfl)%num = iflnum(nfl) ! float identity number cfs(nfl)%type = ifltyp(nfl) ! type of float cfs(nfl)%layer = kfloat(nfl) ! number of model layer cfs(nfl)%pos = flt(nfl,1:3) ! float 3D position cfs(nfl)%starttime = flt(nfl,8) ! initiation time cfs(nfl)%endtime = flt(nfl,9) ! termination time cfs(nfl)%wdepth = flt(nfl,4) ! water-depth cfs(nfl)%temp = flt(nfl,5) ! temperature cfs(nfl)%salt = flt(nfl,6) ! salinity cfs(nfl)%dens = flt(nfl,7) ! density end do call initializecfmod(rt,nflt,rungen,fltflg) print*, 'Calanus individuals initialized before restart' end if call read_daily_carbon(rt) #endif /* CALANUS */ if (mnproc.eq.1) print*, 'Before reading restart, fltflg:', fltflg if (fltflg==2) then ! read restart file - the same as output file c write(lp,*) 'reading restart float parameters' flnmflto = & rungen//'FLT_y'//rt%cyy//'_d'//rt%cdd//'_h'//rt%chh//'.uf' print*, flnmflto open(unit=uoff+801,file=flnmflto,status='unknown', & form='unformatted',action='read') read(uoff+801) new_nflt, itmp, itmp if (mnproc.eq.1) print*, 'new_nflt', new_nflt do nfl=1,new_nflt if (samplfl) then if (ifltyp(nfl).ne.4) then read(uoff+801) iflnum(nfl),rtmp,kfloat(nfl), & (flt(nfl,j),j=1,7) else read(uoff+801) iflnum(nfl),rtmp,kfloat(nfl), & rtmp,rtmp, & rtmp,flt(nfl,10), & (flt(nfl,j),j=4,7), & (flt(nfl,j),j=11,13) endif else if (ifltyp(nfl).ne.4) then read(uoff+801) iflnum(nfl),rtmp,kfloat(nfl), & (flt(nfl,j),j=1,3) else read(uoff+801) iflnum(nfl),rtmp,kfloat(nfl), & rtmp,rtmp, & rtmp,flt(nfl,10) endif endif enddo if (new_nflt < nflt) then do nfl=new_nflt+1, nflt flt(nfl,1)=-999.0 flt(nfl,2)=-999.0 kfloat(nfl)=-1 end do end if close(uoff+801) #if defined (CALANUS) !!! if (mnproc.eq.1) then call read_restart_calanus(rungen,rt,nflt) if (mnproc.eq.1) & call write_calanus_to_file_new(rt,nflt,rungen) c Transfer calanus datato the floats !fix!!!sjekk hva som skjer uten restart file do nfl=1,nflt iflnum(nfl) = cfs(nfl)%num ! float identity number ifltyp(nfl) = cfs(nfl)%type ! type of float kfloat(nfl) = cfs(nfl)%layer ! number of model layer flt(nfl,1:3) = cfs(nfl)%pos ! float 3D position flt(nfl,8) = cfs(nfl)%starttime ! initiation time flt(nfl,9) = cfs(nfl)%endtime ! termination time flt(nfl,4) = cfs(nfl)%wdepth ! water-depth flt(nfl,5) = cfs(nfl)%temp ! temperature flt(nfl,6) = cfs(nfl)%salt ! salinity flt(nfl,7) = cfs(nfl)%dens ! density end do !!! end if #endif /* CALANUS */ end if ! end fltflg==2: read restart file c if (mnproc.eq.1) print*, 'After reading restart, fltflg:', fltflg if (mnproc.eq.1) then write(lp,955) nflt write(lp,956) nfladv write(lp,957) nflsam write(lp,958) min(nflt,10) 955 format(' read initial data for',i6,' floats, time step',i9) 956 format(' float advected every',i7,' baroclinic time steps') 957 format(' float sampled every',i8,' baroclinic time steps'/) 958 format(' input data for first',i3,' floats') do nfl=1,min(nflt,10) write(lp,959) nfl,(flt(nfl,i),i=1,9) 959 format(3x,i2,1p,5e13.5/5x,4e13.5/) enddo call flush(lp) endif !1st tile c c --------------------------------------------------------------------------- c --- allocate arrays for floats c --------------------------------------------------------------------------- c allocate( wvelup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), & wveldn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), & dpdxup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dpdyup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dpdxdn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dpdydn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) allocate( wveli( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy, kdm+1), & uold2( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm), & vold2( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm), & wold2u(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm), & wold2d(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm), & dlondx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dlondy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dlatdx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dlatdy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) if (synflt) then !$OMP PARALLEL DO PRIVATE(j,i,k) !$OMP& SCHEDULE(STATIC,jblk) do j=1-nbdy,jj+nbdy do i=1-nbdy,ii+nbdy do k=1,kk uold2( i,j,k )=huge uold2( i,j,k+kk)=huge vold2( i,j,k )=huge vold2( i,j,k+kk)=huge wold2u(i,j,k )=huge wold2u(i,j,k+kk)=huge wold2d(i,j,k )=huge wold2d(i,j,k+kk)=huge wveli( i,j,k) =huge enddo !k wveli(i,j,kk+1)=huge enddo enddo !$OMP END PARALLEL DO endif !synflt c c --------------------------------------------------------------------------- c --- initialize old horizontal velocities for runga-kutta time interpolation c --- calculate dlondx, dlondy, dlatdx, dlatdy required to convert u, v to c --- longitude/time and latitude/time for float advection c --------------------------------------------------------------------------- c margin = nbdy - 1 c do j=1-margin,jj+margin c do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) do k=1,kk uold2(i,j,k )=u(i,j,k,n)+ubavg(i,j,n) uold2(i,j,k+kk)=u(i,j,k,n)+ubavg(i,j,n) enddo dlondx(i,j)=(plon(i ,j)-plon(i-1,j))/scux(i,j) dlatdx(i,j)=(plat(i ,j)-plat(i-1,j))/scux(i,j) if (plat(i,j)-plat(i-1,j).ne.0.0) then nonlatlon=.true. endif c --- linear correction for curvilinear grid dlon=ulon(i,j)-0.5*(plon(i,j)+plon(i-1,j)) if (dlon.gt.0.0) then dlox=(ulon(i+1,j )-ulon(i ,j ))/scpx(i ,j ) dlax=(ulat(i+1,j )-ulat(i ,j ))/scpx(i ,j ) q_ll= dlon/(plon(i,j)-ulon(i ,j )) dlondx(i,j)=dlondx(i,j)+q_ll*(dlox-dlondx(i,j)) dlatdx(i,j)=dlatdx(i,j)+q_ll*(dlax-dlatdx(i,j)) elseif (dlon.lt.0.0) then dlox=(ulon(i ,j )-ulon(i-1,j ))/scpx(i-1,j ) dlax=(ulat(i ,j )-ulat(i-1,j ))/scpx(i-1,j ) q_ll=-dlon/(ulon(i,j)-plon(i-1,j )) dlondx(i,j)=dlondx(i,j)+q_ll*(dlox-dlondx(i,j)) dlatdx(i,j)=dlatdx(i,j)+q_ll*(dlax-dlatdx(i,j)) endif enddo enddo c do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) do k=1,kk vold2(i,j,k )=v(i,j,k,n)+vbavg(i,j,n) vold2(i,j,k+kk)=v(i,j,k,n)+vbavg(i,j,n) enddo dlondy(i,j)=(plon(i,j )*scpx(i,j )- & plon(i,j-1)*scpx(i,j-1))/scv2(i,j) dlatdy(i,j)=(plat(i,j )*scpx(i,j )- & plat(i,j-1)*scpx(i,j-1))/scv2(i,j) if (plon(i,j)-plon(i,j-1).ne.0.0) then nonlatlon=.true. endif c --- linear correction for curvilinear grid dlat=vlat(i,j)-0.5*(plat(i,j)+plat(i,j-1)) if (dlat.gt.0.0) then dloy=(vlon(i ,j+1)-vlon(i ,j ))/scpy(i ,j ) dlay=(vlat(i ,j+1)-vlat(i ,j ))/scpy(i ,j ) q_ll= dlat/(plat(i,j)-vlat(i ,j )) dlondy(i,j)=dlondy(i,j)+q_ll*(dloy-dlondy(i,j)) dlatdy(i,j)=dlatdy(i,j)+q_ll*(dlay-dlatdy(i,j)) elseif (dlat.lt.0.0) then dloy=(vlon(i ,j )-vlon(i ,j-1))/scpy(i ,j-1) dlay=(vlat(i ,j )-vlat(i ,j-1))/scpy(i ,j-1) q_ll=-dlat/(vlat(i,j)-plat(i ,j-1)) dlondy(i,j)=dlondy(i,j)+q_ll*(dloy-dlondy(i,j)) dlatdy(i,j)=dlatdy(i,j)+q_ll*(dlay-dlatdy(i,j)) endif enddo enddo c do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) do k=1,kk wold2u(i,j,k )=0.0 wold2d(i,j,k )=0.0 wold2u(i,j,k+kk)=0.0 wold2d(i,j,k+kk)=0.0 wveli (i,j,k )=0.0 enddo wveli(i,j,kk+1)=0.0 enddo enddo c enddo !j c c ---------------------------------------- c --- turbulent horizontal velocity option c ---------------------------------------- c c --- initialize turbulence constants c c --- tdecri is in inverse days and dtturb is the turbulent time step in days c c --- dtturb is set to the runga-kutta interpolation time step, which must c --- be within 1 and 2 hours c if (turbvel) then dtturb=deltfl/86400.0 uturb0=sqrt(2.*tbvar*tdecri*dtturb) endif c return end subroutine floats_init subroutine floats_restart c c ----------------------------- c --- output float restart file c ----------------------------- c implicit none c integer i,k,l c if (mnproc.eq.1) then c c --- first remove floats that have run aground or left the domain i=0 do l=1,nflt if (flt(l,1).gt.-999.0) then i=i+1 endif enddo !1 c write(lp,*) 'writing new float input file for restart' call flush(lp) c open(unit=uoff+802,file=flnmfltio, & form='formatted') write(uoff+802,970) i write(uoff+802,970) nstepfl do l=1,nflt if (flt(l,1).gt.-999.0) then if (ifltyp(l).ne.2) then write(uoff+802,971) iflnum(l),ifltyp(l),flt(l,8),flt(l,9), & (flt(l,k),k=1,3) else write(uoff+802,971) iflnum(l),ifltyp(l),flt(l,8),flt(l,9), & (flt(l,k),k=1,2),flt(l,7) endif endif enddo !l close(uoff+802) endif !1st tile 970 format(i9) 971 format(i6,i2,2f9.2,3f15.9) return end subroutine floats_restart subroutine floats(m,n,timefl,ioflag) implicit none c integer, parameter :: nfl_debug = -1 !no debugging * integer, parameter :: nfl_debug = 103 !debug float nfl_debug c include 'common_blocks.h' c integer m,n,ioflag real timefl c c --- local variables integer i,j,k,l real tflt,sflt,thlo,thhi,thflt, & uflt1,uflt2,uflt3,uflt4,uflt, & vflt1,vflt2,vflt3,vflt4,vflt, & wflt1,wflt2,wflt3,wflt4,wflt, & vkflt,tkflt,skflt,xpos0,xpos1,xpos2,xpos3, & ypos0,ypos1,ypos2,ypos3,depflt,plo,phi,q, & depiso,wvhi,wvlo,uvfctr,qisop,ufltm,vfltm, & alomin,alomax,alamin,alamax #if defined (CALANUS) real iflt,max_chl integer max_chl_dep,kcount #endif real p_bottom logical is_inside real uturb,vturb,uturb1,vturb1,dlodx,dlody,dladx,dlady real*4 rtab(200,2) integer ifltll(3),jfltll(3),ngood(3) integer kflt,k1,k2,k3,ier,ngrid,kold1w,kold2w,kold1,kold2, & ntermn,i1,j1,ngoodi integer ied,iede,ifladv c integer*4 seed,numran,inisee,iflag integer*2 iseed1,iseed2 equivalence ( iseed1, iseed2, seed ) c c --- local velocity component storage for synthetic moorings real fltloc(nfldim,3) c c --- 2-d local arrays used for the 16-point box horizontal interpolation real varb2d(4,4),ptlon(4,4,3),ptlat(4,4,3) logical maskpt(4,4,3),maskpi(4,4) #if defined (CALANUS) CAS calanus type (calfin) :: tmpcal real cmirr(nfldim) real cmixl(nfldim) real cphy_conc(nfldim) real cdepth_max_phy_conc(nfldim) real m_nterm(nfldim),m_num(nfldim), ! mortality information & m_day(nfldim) ! to be written to file #endif c c --- arrays required for parallelization !!! real flt1(12*nflt),flt3(17*nflt) real, allocatable :: flt1(:), flt3(:) integer iproc c real timer1, timer2, totime real*8 wtime cAS -- for sorting real tmpflt(13), tmpfltloc(3) integer tmpkfloat, tmpiflnum, tmpifltyp integer new_nflt real, parameter :: dmixl=7.0 real randepth logical, parameter :: diuvm=.false. c include 'stmt_fns.h' c c ------------------------------------------------------------------------- c --- floats.f (hycom version 2.2) c ------------------------------------------------------------------------- c c --- synthetic floats, drifters, and mooring instruments c c --- optionally samples time series of dynamical/thermodynamical variables c --- at the location of each float c c --- four float types are presently supported: c --- 3-d lagrangian (vertically advected by diagnosed vertical velocity) c --- isopycnic (remains on specified density surface) c --- isobaric (remains at the pressure depth where it was released) c --- stationary (synthetic instrument/mooring) c c --- horizontal advection, along with vertical advection of lagrangian floats, c --- is performed using the MICOM runga-kutta-4 time interpolation algorithm c --- developed by Zulema Garraffo c c --- works on any curvilinear grid c c --- float position is stored as longitude and latitude c c --- to horizontal advect the float, u and v are first converted to c --- d(longitude)/dt and d(latitude)/dt as follows: c c --- u_lon = u*dlondx + v*dlondy c --- v_lat = u*dlatdx + v*dlatdy c c --- to horizontally advect the float, terms u*dlondx and u*dlatdx are c --- calculated on u grid points while terms v*dlondy and v*dlatdy are c --- calculated on v grid ponts. these terms are then interpolated to c --- the float locations from the surrounding u and v grid points, c --- respectively c c --- since the terms v*dlondy and u*dlatdx are always zero when the model c --- x,y axes are everywhere lines of constant latitude and longitude, c --- respectively, logical variable 'nonlatlon' prevents horizontal c --- interpolation of these terms in this case c c --- horizontal interpolation to the float location follows the MICOM c --- procedure of Zulema Garraffo - second-order interpolation from the 16 c --- surrounding grid points is performed unless fewer than nptmin (presently c --- set to 10 in subroutine intrph) water points are available, in which c --- case nearest-neighbor interpolation is used. c c --- variables are horizontally interpolated from their native grid (p, u, c --- or v) c c --- adapted for hycom by george halliwell c --- code parallelized by remy baraille c c --- the second-order interpolation subroutine included here is the one c --- included in the mariano and brown parameter matrix objective analysis c --- algorithm to estimate the large scale trend surface - it is different c --- from the algorithm used by Garraffo in MICOM c c --- float initialization is performed in floats_init - information c --- about the input file containing initial float information is c --- presented in that subroutine. c c --- variables in array flt(nfl,n) are: c c --- n = 1 longitude c --- n = 2 latitude c --- n = 3 float depth c --- n = 4 water depth c --- n = 5 temperature c --- n = 6 salinity c --- n = 7 water density c --- n = 8 float start time (in days from start of model run) c --- n = 9 float end time (in days from start of model run) c c --- time series of several variables are interpolated to the location of each c --- float and saved every nflsam time steps when ioflag is set to 1 c c --- variables output onto file float.out for each float are: c c --- 1. initial sequential float number c --- 2. time (in days from start of model run) c --- 3. model layer number that contains the float c --- 4. longitude (u for synthetic moorings) c --- 5. latitude (v for synthetic moorings) c --- 6. float depth (w for synthetic moorings) c --- 7. water depth c --- 8. temperature c --- 9. salinity c --- 10. water density (remains constant for isopycnic floats) c call xctilr(dp( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps) call xctilr(dpu( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vs) call xctilr(dpv( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vs) call xctilr(u( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_uv) call xctilr(v( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vv) call xctilr(ubavg( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_uv) call xctilr(vbavg( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_vv) call xctilr(temp( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps) call xctilr(saln( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps) call xctilr(th3d( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps) c c --- calculate vertical velocity at interfaces as the sum of the c --- vertical interface velocity estimated in cnuity.f and the c --- advective component due to flow parallel to interfaces. c --- wvelup and wveldn represent velocity at the top and bottom of c --- layer k c margin = nbdy c c --- set pressure array at p points do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) do k=1,kk k1=k+1 p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n) enddo enddo enddo enddo c c --- set pressure array at u and v points c --- calculate dpdx, dpdy at interfaces above and below layer k for c --- estimating wvel within layer k c do k=1,kk c margin = nbdy - 1 c do j=1-margin,jj+margin c do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) pu(i,j,k+1)=pu(i,j,k)+dpu(i,j,k,n) if (wvelfl) then dpdxup(i,j)= & (p(i,j,k )*scpy(i,j)-p(i-1,j,k )*scpy(i-1,j)) & /scu2(i,j) dpdxdn(i,j)= & (p(i,j,k+1)*scpy(i,j)-p(i-1,j,k+1)*scpy(i-1,j)) & /scu2(i,j) endif enddo enddo c do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) pv(i,j,k+1)=pv(i,j,k)+dpv(i,j,k,n) if (wvelfl) then dpdyup(i,j)= & (p(i,j,k )*scpx(i,j)-p(i,j-1,k )*scpx(i,j-1)) & /scv2(i,j) dpdydn(i,j)= & (p(i,j,k+1)*scpx(i,j)-p(i,j-1,k+1)*scpx(i,j-1)) & /scv2(i,j) endif enddo enddo c enddo !j c c --- calculate the vertical velocity arrays c margin = nbdy - 2 c do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) if (wvelfl) then wveli(i,j,k+1)=deltfl*wveli(i,j,k+1)/delt1 if (dp(i,j,k,n).gt.tencm) then wvelup(i,j,k)=-0.5*deltfl* & ((u(i ,j ,k,n)+ubavg(i ,j ,n))*dpdxup(i ,j )+ & (u(i+1,j ,k,n)+ubavg(i+1,j ,n))*dpdxup(i+1,j )+ & (v(i ,j ,k,n)+vbavg(i ,j ,n))*dpdyup(i ,j )+ & (v(i ,j+1,k,n)+vbavg(i ,j+1,n))*dpdyup(i ,j+1)) & /onem+wveli(i,j,k ) wveldn(i,j,k)=-0.5*deltfl* & ((u(i ,j ,k,n)+ubavg(i ,j ,n))*dpdxdn(i ,j )+ & (u(i+1,j ,k,n)+ubavg(i+1,j ,n))*dpdxdn(i+1,j )+ & (v(i ,j ,k,n)+vbavg(i ,j ,n))*dpdydn(i ,j )+ & (v(i ,j+1,k,n)+vbavg(i ,j+1,n))*dpdydn(i ,j+1)) & /onem+wveli(i,j,k+1) else wvelup(i,j,k)=0.0 wveldn(i,j,k)=0.0 endif endif enddo enddo enddo !j c enddo !k c * call xctilr(wvelup(1-nbdy,1-nbdy,1),1, kk, 4,4, halo_ps) * call xctilr(wveldn(1-nbdy,1-nbdy,1),1, kk, 4,4, halo_ps) c c --- set old velocity indices used for interpolation of float position c --- perform full float update every two times this subroutine is called; c --- at the intermediate times, just store old velocity fields for the next c --- call c if (mod(nstepfl,nfladv).eq.nfldta) then kold1 =kk kold1w=kk+1 kold2 =0 kold2w=0 go to 10 else kold1 =0 kold1w=0 kold2 =kk kold2w=kk+1 endif c c --- turbulent horizontal velocity option c c --- to customize the parameter settings in blkdat.input that control c --- the calculated turbulent velocity (tbvar, tdecri), refer to c --- Griffa (1996), "Aplications of stochastic particles models to c --- oceanographic problems" in "Stochastic Modelling in Physical c --- Oceanography", pp. 114-140, Adler, Muller, and Rozovoskii, editors c if (turbvel) then c c --- initialize random seeds and create random number table c --- if iflag=1 (iflag=0), then time() is (is not) used to generate seeds c iflag=1 call system_clock(inisee) seed = 414957000-inisee iede = iseed1 ied = iseed2 if (iede .lt. 0) iede = abs(iede) if (ied .lt. 0) ied = abs(ied) c call rantab_ini(rtab,iseed1,iseed2,iflag) endif c c --- get particle locations from processor 1 allocate(flt1(12*nflt)); if (mnproc.eq.1) then do nfl=1,nflt do i=1,9 flt1( i+12*(nfl-1)) = flt(nfl,i) enddo flt1(10+12*(nfl-1)) = kfloat(nfl) flt1(11+12*(nfl-1)) = iflnum(nfl) flt1(12+12*(nfl-1)) = ifltyp(nfl) enddo endif !1st tile call xcastr(flt1(1:12*nflt), 1) if (mnproc.ne.1) then do nfl=1,nflt do i=1,9 flt(nfl,i) = flt1( i+12*(nfl-1)) enddo kfloat(nfl) = flt1(10+12*(nfl-1)) iflnum(nfl) = flt1(11+12*(nfl-1)) ifltyp(nfl) = flt1(12+12*(nfl-1)) enddo endif !1st tile deallocate(flt1) c c ---------------------- c --- begin float loop c ---------------------- c allocate(flt3(17*nflt)) #if defined (CALANUS) m_nterm=huge m_num=huge m_day=huge #endif do nfl=1,nflt c c --- skip if float has previously run aground or exceeded termination time if(flt(nfl,1).le.-999.) then go to 22 endif c c --- ier = error flag for horizontal interpolation c --- ntermn = float termination flag c --- -10 => reached termination time c --- -5 => exited domain c --- >0 => ran aground c --- ifladv = float horizontal advection flag ier=0 ntermn=0 ifladv=1 c c --- suppress float advection if this is the first time step so that c --- the initial position on the output file is exactly the position c --- specified on the input file if (nstepfl.eq.1) then ifladv=0 endif c c --- check if model time has reached float deployment time c --- when it does, again suppress float advection during first pass if (timefl-flt(nfl,8).lt.-0.001) then go to 22 endif if (kfloat(nfl).lt.0) then ifladv=0 endif c c --- check if model time has reached float termination time if (flt(nfl,9).gt.0.0 .and. & timefl-flt(nfl,9).gt. 0.001) then ntermn=-10 endif c c --------------------------------------------------------------------------- c --- search for the processor controlling the tile containing the grid point c --- search for the surrounding 16 grid points on the p, u, and v grids c --------------------------------------------------------------------------- c c --- processeur sur lequel va se derouler le calcul c c AS 08072009 do i=1,3 ifltll(i)=0 jfltll(i)=0 end do c AS 08072009 - end iproc=0 c margin = nbdy - 6 c do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) c c --- search for the surrounding p-grid points if (i+i0.le.itdm .and. j+j0.le.jtdm) then alomin=min(plon(i ,j),plon(i ,j+1), & plon(i+1,j),plon(i+1,j+1)) alomax=max(plon(i ,j),plon(i ,j+1), & plon(i+1,j),plon(i+1,j+1)) alamin=min(plat(i ,j),plat(i ,j+1), & plat(i+1,j),plat(i+1,j+1)) alamax=max(plat(i ,j),plat(i ,j+1), & plat(i+1,j),plat(i+1,j+1)) CAS240709 problem with longitude at the poles: if ((alomax-alomin).gt.100.0) then call correct_longitude(alomin,alomax,flt(nfl,1)) end if CAS240709 CAS: this only work in a regular grid: if (flt(nfl,1).ge.alomin.and.flt(nfl,1).lt.alomax.and. & flt(nfl,2).ge.alamin.and.flt(nfl,2).lt.alamax) then ! check if the point is inside the square call check_square(flt(nfl,1:2), & plon(i:i+1,j:j+1),plat(i:i+1,j:j+1),is_inside) c --- determine if float has exited domain if (i+i0.lt.2 .or. i+i0.gt.itdm-1 .or. & j+j0.lt.2 .or. j+j0.gt.jtdm-1) then ntermn=-5 print*, 'Float number ', nfl , 'will be terminted' endif if (is_inside) then ifltll(1)=i jfltll(1)=j c c --- search for the surrounding u-grid points do i1=i-1,i+1 do j1=j-1,j+1 if (i1+i0.le.itdm .and. j1+j0.le.jtdm) then alomin=min(ulon(i1 ,j1),ulon(i1 ,j1+1), & ulon(i1+1,j1),ulon(i1+1,j1+1)) alomax=max(ulon(i1 ,j1),ulon(i1 ,j1+1), & ulon(i1+1,j1),ulon(i1+1,j1+1)) alamin=min(ulat(i1 ,j1),ulat(i1 ,j1+1), & ulat(i1+1,j1),ulat(i1+1,j1+1)) alamax=max(ulat(i1 ,j1),ulat(i1 ,j1+1), & ulat(i1+1,j1),ulat(i1+1,j1+1)) CAS240709 problem with longitude at the poles: if ((alomax-alomin).gt.100.0) then call correct_longitude(alomin,alomax, & flt(nfl,1)) end if CAS240709 if (flt(nfl,1).ge.alomin .and. & flt(nfl,1).lt.alomax .and. & flt(nfl,2).ge.alamin .and. & flt(nfl,2).lt.alamax) then ifltll(2)=i1 jfltll(2)=j1 call check_square(flt(nfl,1:2), & ulon(i1:i1+1,j1:j1+1),ulat(i1:i1+1,j1:j1+1), & is_inside) c --- determine if float has exited domain if (i1+i0.lt.2 .or. i1+i0.gt.itdm-1 .or. & j1+j0.lt.2 .or. j1+j0.gt.jtdm-1) then ntermn=-5 print*, 'Float number ', nfl , 'will be terminted' go to 11 endif if (is_inside) then ifltll(2)=i1 jfltll(2)=j1 go to 11 endif !is_inside endif !alomin/max endif !itdm/jtdm enddo !j1 enddo !i1 11 continue c c --- search for the surrounding v-grid points do i1=i-1,i+1 do j1=j-1,j+1 if (i1+i0.le.itdm .and. j1+j0.le.jtdm) then alomin=min(vlon(i1 ,j1),vlon(i1 ,j1+1), & vlon(i1+1,j1),vlon(i1+1,j1+1)) alomax=max(vlon(i1 ,j1),vlon(i1 ,j1+1), & vlon(i1+1,j1),vlon(i1+1,j1+1)) alamin=min(vlat(i1 ,j1),vlat(i1 ,j1+1), & vlat(i1+1,j1),vlat(i1+1,j1+1)) alamax=max(vlat(i1 ,j1),vlat(i1 ,j1+1), & vlat(i1+1,j1),vlat(i1+1,j1+1)) CAS240709 problem with longitude at the poles: if ((alomax-alomin).gt.100.0) then call correct_longitude(alomin,alomax, & flt(nfl,1)) end if CAS240709 if (flt(nfl,1).ge.alomin .and. & flt(nfl,1).lt.alomax .and. & flt(nfl,2).ge.alamin .and. & flt(nfl,2).lt.alamax) then call check_square(flt(nfl,1:2), & vlon(i1:i1+1,j1:j1+1),vlat(i1:i1+1,j1:j1+1), & is_inside) c --- determine if float has exited domain if (i1+i0.lt.2 .or. i1+i0.gt.itdm-1 .or. & j1+j0.lt.2 .or. j1+j0.gt.jtdm-1) then print*, 'Float number ', nfl , 'will be terminted' ntermn=-5 go to 12 endif if (is_inside) then ifltll(3)=i1 jfltll(3)=j1 go to 12 endif !is_inside endif !alomin/max endif !itdm/jtdm enddo !j1 enddo !i1 12 continue c c --- set processor number for the tile containing the float and exit grid c --- point search iproc=mnproc go to 13 endif !is_inside endif !alomin/max endif !itdm/jtdm enddo !i enddo !l enddo !j c 13 continue c c --- float nfl is now updated by the processor running the tile containing c --- the float c if (iproc.gt.0) then c if (nfl.eq.nfl_debug) then write(lp,100) nfl,nflt,nstep,nstepfl write(lp,101) nfl,ntermn,flt(nfl,1),flt(nfl,2),flt(nfl,3), & (ifltll(i)+i0,jfltll(i)+j0,i=1,3), & mnproc, & plon(ifltll(1),jfltll(1)), & plat(ifltll(1),jfltll(1)) 100 format(/'diagnostics for float',i6,' of',i6,', time step', & i9/'float time step',i9) 101 format('float',i6,' ntermn:',i6,' position:',3(1pe12.4)/ & 'lower left points, p,u,v:',3(2i5,2x)/ & 'mnproc =',i4,' plon,plat =',1pe12.4,1pe12.4) call flush(lp) endif !nfl_debug c c --- if float has exited domain or run aground as determined previously, c --- jump ahead to the float termination code bloci if (ntermn.eq.-5 .or. ntermn.eq.1) then go to 30 endif c c --- set ptlon, ptlat for all horizontal interpolations from each grid ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 ptlon(i1,j1,ngrid)=plon(i,j) ptlat(i1,j1,ngrid)=plat(i,j) enddo enddo c ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 ptlon(i1,j1,ngrid)=ulon(i,j) ptlat(i1,j1,ngrid)=ulat(i,j) enddo enddo c ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 ptlon(i1,j1,ngrid)=vlon(i,j) ptlat(i1,j1,ngrid)=vlat(i,j) enddo enddo c c --- the float is assumed to remain within the same layer that it was in c --- during the previous update unless it is the first advection time step c --- for the float, in which case it is initially assumed to be in layer 1 c k=max(1,kfloat(nfl)) c c --- mask p grid points that are on land or where the layer containing the c --- float is a zero or nearly-zero thickness layer at the bottom ngrid=1 ngood(ngrid)=0 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (depths(i,j).lt.0.01 .or. depths(i,j).eq.huge .or. & depths(i,j)*onem-p(i,j,k).lt.tencm) then maskpt(i1,j1,ngrid)=.true. else maskpt(i1,j1,ngrid)=.false. ngood(ngrid)=ngood(ngrid)+1 endif enddo enddo c if (nfl.eq.nfl_debug) then write(lp,102) ngood(1) 102 format('initial masking: no. of good p points',i4) endif !nfl_debug c c -------------------------------------------------- c --- determine the model layer containing the float c -------------------------------------------------- c c --- if this is not the first advection time step for the float, determine c --- if the float has transited to another model layer c if (kfloat(nfl).ge.1) then c c --- check if the float is deeper than the interpolated bottom c --- depth - if so, reset the depth to 0.1 m above the bottom - this c --- prevents from running aground too frequently in coastal regions if (ifltyp(nfl).ne.4) then ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=p(i,j,kk+1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,p_bottom,ier) flt(nfl,3)=max(fldepm, & min(flt(nfl,3),(p_bottom-onem)/onem)) endif c c --- has the float remained in the same layer? ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=p(i,j,k) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,phi,ier) c c --- terminate float because it has run aground if (ier.eq.999) then ntermn=2 go to 30 endif c do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=p(i,j,k+1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,plo,ier) phi=phi/onem plo=plo/onem c c --- if the float is in the same model layer, skip ahead if (flt(nfl,3).gt.phi .and. flt(nfl,3).le.plo) then kflt=kfloat(nfl) go to 40 endif c c --- determine the new model layer containing the float c c --- did the float move up one layer? if (flt(nfl,3).le.phi .and. kfloat(nfl).gt.1) then k=kfloat(nfl)-1 plo=phi do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=p(i,j,k) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,phi,ier) phi=phi/onem if (flt(nfl,3).gt.phi .and. flt(nfl,3).le.plo) then kflt=k go to 40 endif endif c c --- did the float move down one layer? if (flt(nfl,3).gt.plo .and. kfloat(nfl).lt.kk) then k=kfloat(nfl)+1 phi=plo do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=p(i,j,k) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,plo,ier) plo=plo/onem if (flt(nfl,3).gt.phi .and. flt(nfl,3).le.plo) then kflt=k go to 40 endif endif c endif c c --- if the float did not move up or down one layer, or if this is the c --- first time step after float initialization, search for the layer c --- containg the float from top to bottom ngrid=1 plo=0. k2=0 k3=0 do k=1,kk k1=min(k+3,kk+1) c c --- if the depth of interface k1 is not deeper than the float at any c --- surrounding grid points, skip to end of k-loop if (k2.eq.0) then do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (.not.maskpt(i1,j1,ngrid) .and. & p(i,j,k1)/onem.gt.flt(nfl,3)) then k2=1 endif enddo enddo endif if (k2.eq.0) then go to 42 endif c c --- mask points where the layer being tested rests on the bottom do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (k2.eq.1) then maskpi(i1,j1)=maskpt(i1,j1,ngrid) endif if (.not.maskpi(i1,j1) .and. & depths(i,j)*onem-p(i,j,k).lt.tencm) then maskpi(i1,j1)=.true. ngood(ngrid)=ngood(ngrid)-1 endif enddo enddo k2=2 c c --- if too few good points remain, float has run aground if (ngood(ngrid).le.2) then ntermn=3 go to 30 endif c c --- interpolate to find plo phi=plo ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=p(i,j,k+1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,plo,ier) plo=plo/onem if(plo-phi.gt.0.1) then k3=k else kflt=k3 do i1=1,4 do j1=1,4 maskpt(i1,j1,ngrid)=maskpi(i1,j1) enddo enddo go to 40 endif if (flt(nfl,3).lt.plo) then kflt=k do i1=1,4 do j1=1,4 maskpt(i1,j1,ngrid)=maskpi(i1,j1) enddo enddo go to 40 endif 42 continue enddo c c --- if layer selection fails, assume float has run aground print *,'warning - selection of float layer failed' ntermn=4 go to 30 c c --- we now know the model layer containing the float 40 kfloat(nfl)=kflt c c ------------------------------------------ c --- set grid masks for u and v grid points c ------------------------------------------ c k=kflt ngrid=2 ngood(ngrid)=0 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (depthu(i,j).lt.onecm .or. depthu(i,j).eq.huge .or. & depthu(i,j)-pu(i,j,k).lt.tencm) then maskpt(i1,j1,ngrid)=.true. else maskpt(i1,j1,ngrid)=.false. ngood(ngrid)=ngood(ngrid)+1 endif enddo enddo c ngrid=3 ngood(ngrid)=0 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (depthv(i,j).lt.onecm .or. depthv(i,j).eq.huge .or. & depthv(i,j)-pv(i,j,k).lt.tencm) then maskpt(i1,j1,ngrid)=.true. else maskpt(i1,j1,ngrid)=.false. ngood(ngrid)=ngood(ngrid)+1 endif enddo enddo c c --- if too few good p, u, v grid points, assume float has run aground if (ngood(1).le.2 .or. ngood(2).le.2 .or. ngood(3).le.2) then ntermn=5 go to 30 endif c if (nfl.eq.nfl_debug) then write(lp,103) nfl,k,phi,flt(nfl,3),plo,q write(lp,104) ngood 103 format('nfl,k,phi,pflt,plo,q',i6,i3,1p,4e12.4) 104 format('number of good points, p,u,v',3i4) endif !nfl_debug c c --- set vertical indices, and also set q for vertical interpolation k=kfloat(nfl) k1=max(1,k-1) k2=min(kk,k+1) q=(plo-flt(nfl,3))/max(0.001,plo-phi) q=max(q,0.0) ! AS q has som outragous values q=min(q,1.0) c c ---------------------------------------------------------- c --- advect the float horizontally, then move it vertically c --- do not move float if it is a mooring instrument c ---------------------------------------------------------- c c ---------------------------------- c --- horizontal advection of floats c ---------------------------------- c if (ifltyp(nfl).ne.4) then c c --- interpolate u,v to float location - velocities at the new time and c --- at two earlier times are used for the runga-kutta time interpolation c xpos0=flt(nfl,1) ypos0=flt(nfl,2) c ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondx(i,j)*uold2(i,j,k+kold2) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,uflt1,ier) c c --- terminate float because it has run aground if (ier.eq.999) then ntermn=6 go to 30 endif c if (nonlatlon) then ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondy(i,j)*vold2(i,j,k+kold2) enddo enddo if (l.eq.1) then call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,ufltm,ier) uflt1=uflt1+ufltm endif c c --- terminate float because it has run aground if (ier.eq.999) then ntermn=7 go to 30 endif endif c ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdy(i,j)*vold2(i,j,k+kold2) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vflt1,ier) c c --- terminate float because it has run aground if (ier.eq.999) then ntermn=8 go to 30 endif c if (nonlatlon) then ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdx(i,j)*uold2(i,j,k+kold2) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vfltm,ier) vflt1=vflt1+vfltm c c --- terminate float because it has run aground if (ier.eq.999) then ntermn=9 go to 30 endif endif c xpos1=flt(nfl,1)+uflt1 ypos1=flt(nfl,2)+vflt1 c ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondx(i,j)*uold2(i,j,k+kold1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos1,ypos1,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,uflt2,ier) c if (nonlatlon) then ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondy(i,j)*vold2(i,j,k+kold1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,ufltm,ier) uflt2=uflt2+ufltm endif c ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdy(i,j)*vold2(i,j,k+kold1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos1,ypos1,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vflt2,ier) c if (nonlatlon) then ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdx(i,j)*uold2(i,j,k+kold1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vfltm,ier) vflt2=vflt2+vfltm endif c xpos2=flt(nfl,1)+uflt2 ypos2=flt(nfl,2)+vflt2 ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondx(i,j)*uold2(i,j,k+kold1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos2,ypos2,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,uflt3,ier) c if (nonlatlon) then ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondy(i,j)*vold2(i,j,k+kold1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,ufltm,ier) uflt3=uflt3+ufltm endif c ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdy(i,j)*vold2(i,j,k+kold1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos2,ypos2,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vflt3,ier) c if (nonlatlon) then ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdx(i,j)*uold2(i,j,k+kold1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vfltm,ier) vflt3=vflt3+vfltm endif c endif c if (ifltyp(nfl).ne.4) then xpos3=flt(nfl,1)+2.0*uflt3 ypos3=flt(nfl,2)+2.0*vflt3 else xpos3=flt(nfl,1) ypos3=flt(nfl,2) endif c ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (ifltyp(nfl).ne.4) then varb2d(i1,j1)=dlondx(i,j)*(u(i,j,k,n)+ubavg(i,j,n)) else varb2d(i1,j1)=u(i,j,k,n)+ubavg(i,j,n) endif enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos3,ypos3,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,uflt4,ier) c if (nonlatlon .and. ifltyp(nfl).ne.4) then ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondy(i,j)*(v(i,j,k,n)+vbavg(i,j,n)) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,ufltm,ier) uflt4=uflt4+ufltm endif c ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (ifltyp(nfl).ne.4) then varb2d(i1,j1)=dlatdy(i,j)*(v(i,j,k,n)+vbavg(i,j,n)) else varb2d(i1,j1)=v(i,j,k,n)+vbavg(i,j,n) endif enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos3,ypos3,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vflt4,ier) c if (nonlatlon .and. ifltyp(nfl).ne.4) then ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdx(i,j)*(u(i,j,k,n)+ubavg(i,j,n)) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vfltm,ier) vflt4=vflt4+vfltm endif c if (ifltyp(nfl).ne.4) then uflt=(uflt1+2.0*uflt2+2.0*uflt3+uflt4)/6.0 vflt=(vflt1+2.0*vflt2+2.0*vflt3+vflt4)/6.0 c c --- add turbulent velocity to u and v if requested c if (turbvel) then call rantab(rtab,iseed1,iseed2,numran) uturb=uturb0*rtab(numran, 2) call rantab(rtab,iseed1,iseed2,numran) vturb=uturb0*rtab(numran, 2) c c --- convert uturb, vturb to d(longitude)/dt, d(latitude)/dt c --- interpolate dlondx, dlondy, dlatdx, dlatdy to the float location c ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondx(i,j) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,dlodx,ier) if (nonlatlon) then do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdx(i,j) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,dladx,ier) ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlondy(i,j) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,dlody,ier) else dladx=0.0 dlody=0.0 endif ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dlatdy(i,j) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,dlady,ier) c c --- add turbulent velocity to uflt, vflt c uturb1=uturb*dlodx+vturb*dlody vturb1=uturb*dladx+vturb*dlady c if (nfl.eq.nfl_debug) then write(lp,105) nfl,uflt,uturb1,vflt,vturb1 105 format('nfl,uflt,uturb1,vflt,vturb1',i5,1p,4e12.4) endif !nfl_debug c uflt=(1.0-dtturb*tdecri)*uflt+uturb1 vflt=(1.0-dtturb*tdecri)*vflt+vturb1 c endif c c --- update float horizontal position c if (ifladv.eq.1) then flt(nfl,1)=flt(nfl,1)+2.0*deltfl*uflt flt(nfl,2)=flt(nfl,2)+2.0*deltfl*vflt endif c if (nfl.eq.nfl_debug) then write(lp,106) nfl,flt(nfl,1),xpos1,xpos2,xpos3, & uflt1,uflt2,uflt3,uflt4,2.0*deltfl*uflt write(lp,107) nfl,flt(nfl,2),ypos1,ypos2,ypos3, & vflt1,vflt2,vflt3,vflt4,2.0*deltfl*vflt write(lp,108) nfl,uflt,vflt,flt(nfl,1),flt(nfl,2) 106 format('nfl,x-position',i6,1p,4e12.4/' u-velocity',6x,5e12.4) 107 format('nfl,y-position',i6,1p,4e12.4/' v-velocity',6x,5e12.4) 108 format('nfl,udeg,vdeg',i6,1p,2e13.5/'new position',2e13.5) endif !nfl_debug c endif c c ----------------------------------------------- c --- vertical advection of 3-d lagrangian floats c ----------------------------------------------- c c --- horizontally interpolate diagnosed vertical velocity to the float c --- location at the new time and at two earlier times to execute the c --- runga-kutta time interpolation c if (wvelfl .and. ifltyp(nfl).eq.1) then c ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 wvhi=wold2u(i,j,k+kold2) wvlo=wold2d(i,j,k+kold2) varb2d(i1,j1)=q*wvhi+(1.0-q)*wvlo enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos0,ypos0,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,wflt1,ier) c ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 wvhi=wold2u(i,j,k+kold1) wvlo=wold2d(i,j,k+kold1) varb2d(i1,j1)=q*wvhi+(1.0-q)*wvlo enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos1,ypos1,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,wflt2,ier) c call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos2,ypos2,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,wflt3,ier) c endif c if (wvelfl) then ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 wvhi=wvelup(i,j,k) wvlo=wveldn(i,j,k) varb2d(i1,j1)=q*wvhi+(1.0-q)*wvlo enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & xpos3,ypos3,maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,wflt4,ier) !ASD if (wflt4>200) then write(*,'(a,i5)') 'Float ', nfl write(*,'(a,f9.3)') ' wflt: ', wflt4 write(*,'(a,f9.3)') 'Water depth: ', flt(nfl,4) write(*,'(a,4f11.4,i5)') 'Var2d:', varb2d(:,1) write(*,'(a,4f11.4,i5)') 'Var2d:', varb2d(:,2) write(*,'(a,4f11.4,i5)') 'Var2d:', varb2d(:,3) write(*,'(a,4f11.4,i5)') 'Var2d:', varb2d(:,4) write(*,'(a,4f8.4,i5)') 'ptlon:', ptlon(:,1,ngrid),nfl write(*,'(a,4f8.4,i5)') 'ptlon:', ptlon(:,2,ngrid),nfl write(*,'(a,4f8.4,i5)') 'ptlon:', ptlon(:,3,ngrid),nfl write(*,'(a,4f8.4,i5)') 'ptlon:', ptlon(:,4,ngrid),nfl write(*,'(a,4f8.4,i5)') 'ptlat:', ptlat(:,1,ngrid),nfl write(*,'(a,4f8.4,i5)') 'ptlat:', ptlat(:,2,ngrid),nfl write(*,'(a,4f8.4,i5)') 'ptlat:', ptlat(:,3,ngrid),nfl write(*,'(a,4f8.4,i5)') 'ptlat:', ptlat(:,4,ngrid),nfl write(*,'(a,l4,i5)') 'maskpt:', maskpt(:,1,ngrid),nfl write(*,'(a,l4,i5)') 'maskpt:', maskpt(:,2,ngrid),nfl write(*,'(a,l4,i5)') 'maskpt:', maskpt(:,3,ngrid),nfl write(*,'(a,l4,i5)') 'maskpt:', maskpt(:,4,ngrid),nfl write(*,'(a,f8.4)') 'xpos3: ', xpos3 write(*,'(a,f8.4)') 'ypos3: ', ypos3 write(*,'(a,i5)') 'ngood: ', ngood(ngrid) write(*,'(a,i5)') 'ngrid: ', ngrid write(*,'(a,i5)') 'intpfl: ', intpfl write(*,'(a,f8.4)') 'radian: ', radian write(*,'(a,f9.4)') 'wflt4: ', wflt4 write(*,'(a,i5)') 'ier: ', ier end if !ASD endif c if (ifltyp(nfl).eq.1) then c c --- dividing by 3 gives total vertical displacement over interval 2*deltfl wflt=(wflt1+2.0*wflt2+2.0*wflt3+wflt4)/3.0 if (ifladv.eq.1) then !AS 07292009: This may cause floats to aggregate at the surface ! flt(nfl,3)=max(fldepm,flt(nfl,3)-wflt) !AS 07292009: Try reflective boundary condition: flt(nfl,3)=flt(nfl,3)-wflt if (flt(nfl,3)<=dmixl) then call random_number(randepth) flt(nfl,3) = min(dmixl*randepth,(p_bottom-onem)/onem) endif endif if (nfl.eq.nfl_debug) then write(lp,109) nfl,wflt,flt(nfl,3) 109 format('nfl,w,new depth',i6,1p,2e13.5) endif !nfl_debug c elseif (ifltyp(nfl).eq.2) then c c --------------------------------------- c --- vertical motion of isopycnic floats c --------------------------------------- c c --- set float depth to the vertically interpolated depth of the specified c --- density interface c c --- assume float runs aground if no water denser than the specified density c --- exists in the water column c c --- limit float depth to fldepm if no water lighter than the specified c --- density exists in the water column c c --- find smallest k where float density is exceeded by the density in c --- layer k at at least one of the good 16 surrounding grid points thflt=flt(nfl,7)-thbase ngrid=1 ngoodi=ngood(ngrid) do k=1,kk do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 maskpi(i1,j1)=maskpt(i1,j1,ngrid) if (.not.maskpi(i1,j1)) then if (th3d(i,j,k,n).gt.thflt) then kflt=k go to 33 endif endif enddo enddo enddo 33 continue c c --- move downward from layer kflt to find the two layers with densities c --- bracketing the assigned float density do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (.not.maskpi(i1,j1)) then if (kflt.eq.1) then phi=0. plo=.5*dp(i,j,kflt,n)/onem thflt=flt(nfl,7)-thbase if(thflt.le.th3d(i,j,kflt,n)) then varb2d(i1,j1)=max(fldepm,plo) go to 50 endif endif do k=max(2,kflt),kk phi=0.5*(p(i,j,k)+p(i,j,k-1))/onem plo=0.5*(p(i,j,k+1)+p(i,j,k))/onem thhi=th3d(i,j,k-1,n) thlo=th3d(i,j,k ,n) if(thflt.gt.thhi .and. thflt.le.thlo) then qisop=max(0.0,min(1.0,(thlo-thflt)/ & max(epsil,thlo-thhi))) if (nfl.eq.nfl_debug) then write(lp,*) 'QISOP', thlo, thflt, thhi, qisop write(lp,*) 'DEPTH', plo, qisop*phi+(1.-qisop)*plo, phi end if varb2d(i1,j1)=max(fldepm,qisop*phi+(1.-qisop)*plo) if (depths(i,j)-varb2d(i1,j1).lt.0.1) then maskpi(i1,j1)=.true. ngoodi=ngoodi-1 endif go to 50 endif enddo 50 continue endif enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpi, & ngoodi,ngrid,intpfl,radian,depiso,ier) flt(nfl,3)=max(fldepm,depiso) if (nfl.eq.nfl_debug) then write(lp,110) nfl,flt(nfl,3), max(fldepm,depiso) 110 format('nfl,new isopycnic depth',i6,1p,2e13.5) endif !nfl_debug c elseif (ifltyp(nfl).eq.4) then c wflt=wflt4/deltfl elseif (diuvm.eqv..true.) then if (flt(nfl,3).ne.maxdepth(nfl)) then if (rt%ihh.ge.6.and.rt%ihh.lt.18) then call dvm(flt(nfl,3),baclin*4.0,0, & maxdepth(nfl),mindepth(nfl)) endif endif if (flt(nfl,3).ne.mindepth(nfl)) then if (rt%ihh.lt.6.or.rt%ihh.ge.18) then call dvm(flt(nfl,3),baclin*4.0,1, & maxdepth(nfl),mindepth(nfl)) endif endif c endif c c -------------------------------------------------------------- c --- interpolate water properties to the present float location c -------------------------------------------------------------- c if (ioflag.eq.1 .and. samplfl) then c c --- recalculate q. assume float does not leave the layer that it was c --- originally diagnosed to be within. c q=max(0.0,min(1.0,(plo-flt(nfl,3))/max(0.001,plo-phi))) c c --- water depth interpolation if (ifltyp(nfl).ne.4 .or. flt(nfl,4).eq.0.0) then ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=depths(i,j) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & 16,ngrid,intpfl,radian,depflt,ier) flt(nfl,4)=depflt endif c c --- temperature interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (ifltyp(nfl).eq.4) then k1=max(1 ,k-1) k2=min(kk,k+1) varb2d(i1,j1)= & 0.5*( q *(temp(i,j,k1,n)+temp(i,j,k,n))+ & (1.0-q)*(temp(i,j,k2,n)+temp(i,j,k,n))) else varb2d(i1,j1)=temp(i,j,k,n) endif enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,tflt,ier) flt(nfl,5)=tflt c c --- salinity interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (ifltyp(nfl).eq.4) then k1=max(1 ,k-1) k2=min(kk,k+1) varb2d(i1,j1)= & 0.5*( q *(saln(i,j,k1,n)+saln(i,j,k,n))+ & (1.0-q)*(saln(i,j,k2,n)+saln(i,j,k,n))) else varb2d(i1,j1)=saln(i,j,k,n) endif enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,sflt,ier) flt(nfl,6)=sflt c c --- calculate density if (ifltyp(nfl).ne.2) then thflt=sig(tflt,sflt) flt(nfl,7)=thflt else tflt=tofsig(flt(nfl,7),sflt) flt(nfl,5)=tflt endif #if defined (CALANUS) C --- irradiance interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=radmax(i,j) enddo enddo call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), & ngood(ngrid),ngrid,intpfl,radian,iflt,ier) cmirr(nfl)=iflt C --- mixed layer interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dpmixl(i,j,n)/onem enddo enddo call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), & ngood(ngrid),ngrid,intpfl,radian,iflt,ier) cmixl(nfl)=iflt C --- phytoplankton carbon interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (ifltyp(nfl).eq.4) then k1=max(1 ,k-1) k2=min(kk,k+1) varb2d(i1,j1)= & 0.5*( q *(mod_carb(i,j,k1)+mod_carb(i,j,k))+ & (1.0-q)*(mod_carb(i,j,k2)+mod_carb(i,j,k))) else varb2d(i1,j1)=mod_carb(i,j,k) endif enddo enddo call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), & ngood(ngrid),ngrid,intpfl,radian,iflt,ier) cphy_conc(nfl)=iflt C --- find depth of maximum phytoplankton ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 max_chl=-1.0 do kcount=1,kdm if (mod_carb(i,j,kcount).gt.max_chl) then max_chl=mod_carb(i,j,kcount) max_chl_dep=kcount end if end do varb2d(i1,j1)=0.5*(p(i,j,max_chl_dep)+ & p(i,j,max_chl_dep+1))/onem enddo enddo call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), & ngood(ngrid),ngrid,intpfl,radian,iflt,ier) cdepth_max_phy_conc(nfl)=iflt #endif c c --- interpolate model fields for float type 4 (stationary) if (ifltyp(nfl).eq.4) then c c --- 3-d u and v interpolation ngrid=2 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 k1=max(1 ,k-1) k2=min(kk,k+1) varb2d(i1,j1)=ubavg(i,j,n)+ & 0.5*( q *(u(i,j,k1,n)+u(i,j,k,n))+ & (1.0-q)*(u(i,j,k2,n)+u(i,j,k,n))) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,uflt,ier) c ngrid=3 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 k1=max(1 ,k-1) k2=min(kk,k+1) varb2d(i1,j1)=vbavg(i,j,n)+ & 0.5*( q *(v(i,j,k1,n)+v(i,j,k,n))+ & (1.0-q)*(v(i,j,k2,n)+v(i,j,k,n))) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vflt,ier) c c --- wveli interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=q*wveli(i,j,k)+(1.0-q)*wveli(i,j,k+1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,wflt1,ier) flt(nfl,10)=wflt1/deltfl c c --- vertical viscosity coefficient interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=q*vcty(i,j,k)+(1.0-q)*vcty(i,j,k+1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,vkflt,ier) flt(nfl,11)=max(1.0e-4,vkflt) c c --- vertical temperature diffusivity coefficient interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=q*dift(i,j,k)+(1.0-q)*dift(i,j,k+1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,tkflt,ier) flt(nfl,12)=max(1.0e-5,tkflt) c c --- vertical scalar diffusivity coefficient interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=q*difs(i,j,k)+(1.0-q)*difs(i,j,k+1) enddo enddo call intrph(varb2d,ptlon(:,:,ngrid),ptlat(:,:,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(:,:,ngrid), & ngood(ngrid),ngrid,intpfl,radian,skflt,ier) flt(nfl,13)=max(1.0e-5,skflt) c endif ! float type 4 c if (nfl.eq.nfl_debug) then write(lp,111) nfl,tflt,sflt,thflt 111 format('new t,s,th, float',i5,2x,3f9.3) endif !nfl_debug c endif ! ioflag.eq.1 #if defined (CALANUS) C --- find properties important for the calanus C --- irradiance interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=radmax(i,j) enddo enddo call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), & ngood(ngrid),ngrid,intpfl,radian,iflt,ier) cmirr(nfl)=iflt C --- mixed layer interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 varb2d(i1,j1)=dpmixl(i,j,n)/onem enddo enddo call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), & ngood(ngrid),ngrid,intpfl,radian,iflt,ier) cmixl(nfl)=iflt C --- phytoplankton carbon interpolation ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 if (ifltyp(nfl).eq.4) then k1=max(1 ,k-1) k2=min(kk,k+1) varb2d(i1,j1)= & 0.5*( q *(mod_carb(i,j,k1)+mod_carb(i,j,k))+ & (1.0-q)*(mod_carb(i,j,k2)+mod_carb(i,j,k))) else varb2d(i1,j1)=mod_carb(i,j,k) endif enddo enddo call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), & ngood(ngrid),ngrid,intpfl,radian,iflt,ier) cphy_conc(nfl)=iflt C --- find depth of maximum phytoplankton ngrid=1 do i1=1,4 i=i1+ifltll(ngrid)-2 do j1=1,4 j=j1+jfltll(ngrid)-2 max_chl=-1.0 do kcount=1,kdm if (mod_carb(i,j,kcount).gt.max_chl) then max_chl=mod_carb(i,j,kcount) max_chl_dep=kcount end if end do varb2d(i1,j1)=0.5*(p(i,j,max_chl_dep)+ & p(i,j,max_chl_dep+1))/onem enddo enddo call intrph(varb2d,ptlon(1,1,ngrid),ptlat(1,1,ngrid), & flt(nfl,1),flt(nfl,2),maskpt(1,1,ngrid), & ngood(ngrid),ngrid,intpfl,radian,iflt,ier) cdepth_max_phy_conc(nfl)=iflt #endif c c ---------------------------------------------------- c --- flag floats that have run aground or left domain c ---------------------------------------------------- c c --- jump ahead if the float has not reached its termination time if (ntermn.gt.-10) go to 20 c c --- the following code block is entered when the previous code has c --- detected that the float needs to be terminated 30 flt(nfl,1)=-999. flt(nfl,2)=-999. if (ntermn.gt.0) then write(lp,112) nfl,nstep,nstepfl,ier,ntermn,ngood 112 format('float',i6,' aground, time step',i9, & ' float time step',i9/' ier ',i3,' exit point', & i3,' ngood 1-4',4i3) #if defined (CALANUS) m_nterm(nfl)=real(ntermn) m_num(nfl)=real(iflnum(nfl)) m_day(nfl)=real(rt%idd) #endif elseif (ntermn.eq.-5) then write(lp,113) nfl,iflnum(nfl) 113 format('float',2i6,' exits domain') #if defined (CALANUS) m_nterm(nfl)=real(ntermn) m_num(nfl)=real(iflnum(nfl)) m_day(nfl)=real(rt%idd) #endif elseif (ntermn.eq.-10) then write(lp,114) nfl 114 format('float',i6,' reaches termination time') #if defined (CALANUS) m_nterm(nfl)=real(ntermn) m_num(nfl)=real(iflnum(nfl)) m_day(nfl)=real(rt%idd) #endif endif c 20 continue c c ---------------------------------------------------- c --- store velocity components for synthetic moorings c ---------------------------------------------------- c if (ifltyp(nfl).eq.4) then fltloc(nfl,1)=uflt fltloc(nfl,2)=vflt fltloc(nfl,3)=wflt else fltloc(nfl,1)=0.0 fltloc(nfl,2)=0.0 fltloc(nfl,3)=0.0 endif c c--- store updated particle location c do i=1,13 flt3( i+17*(nfl-1)) = flt( nfl,i) enddo do i=1,3 flt3(13+i+17*(nfl-1)) = fltloc(nfl,i) enddo flt3(17+17*(nfl-1)) = kfloat(nfl) c c --- else du premier if de la boucle else !float is not on this tile, so disable flt3 do i=1,17 flt3(i+17*(nfl-1)) = huge flt(nfl,1)=-999.0 flt(nfl,2)=-999.0 #if defined (CALANUS) cmirr(nfl) = huge cmixl(nfl) = huge cphy_conc(nfl) = huge cdepth_max_phy_conc(nfl) = huge #endif enddo endif c go to 222 c 22 continue c c --- float is unchanged do i=1,17 flt3(i+17*(nfl-1)) = huge enddo c 222 continue c c --------------------- c --- end of float loop c --------------------- c enddo !nfl c --- copy floats back to 1st processor call xcminr(flt3(1:17*nflt), 1) #if defined (CALANUS) call xcminr(cmirr(1:nflt),1) call xcminr(cmixl(1:nflt),1) call xcminr(cphy_conc(1:nflt),1) call xcminr(cdepth_max_phy_conc(1:nflt),1) call xcminr(m_nterm(1:nflt),1) call xcminr(m_num(1:nflt),1) call xcminr(m_day(1:nflt),1) ! file for mortality reasons if (mnproc.eq.1) then open(unit=811,file='mortality_out',status='old', & form='formatted',position='append') 119 format(a12,i12,a8,i8,a45) do nfl=1,nflt if (int(m_nterm(nfl)).gt.0) write(811,119) 'Individual #', & int(m_num(nfl)), ' on day ', int(m_day(nfl)), & ' killed - float went aground' if (int(m_nterm(nfl)).eq.-5) write(811,119) 'Individual #', & int(m_num(nfl)), ' on day ', int(m_day(nfl)), & ' killed - float exited domain' if (int(m_nterm(nfl)).eq.-10) write(811,119) 'Individual #', & int(m_num(nfl)), ' on day ',int( m_day(nfl)), & ' killed - float reached termination time' end do close(811) endif #endif if (mnproc.eq.1) then do nfl=1,nflt if (flt3(1+17*(nfl-1)).ne.huge) then do i=1,13 flt( nfl,i) = flt3( i+17*(nfl-1)) enddo do i=1,3 fltloc(nfl,i) = flt3(13+i+17*(nfl-1)) enddo kfloat(nfl) = flt3(17+17*(nfl-1)) endif !updated float enddo !nfl endif !1st tile deallocate(flt3) #if defined (CALANUS) CAS update the calanus positions !READ IN DAILY CARBON FROM NORWECOM ONCE A DAY if (.not.idealized) then if (rt%ihh==0.and.rt%iss==0) then call read_daily_carbon(rt) end if end if if (mnproc.eq.1) then do nfl=1,nflt cfs(nfl)%num = iflnum(nfl) ! float identity number cfs(nfl)%type = ifltyp(nfl) ! type of float cfs(nfl)%layer = kfloat(nfl) ! number of model layer cfs(nfl)%pos = flt(nfl,1:3) ! float 3D position cfs(nfl)%starttime = flt(nfl,8) ! initiation time cfs(nfl)%endtime = flt(nfl,9) ! termination time cfs(nfl)%wdepth = flt(nfl,4) ! water-depth cfs(nfl)%temp = flt(nfl,5) ! temperature cfs(nfl)%salt = flt(nfl,6) ! salinity cfs(nfl)%dens = flt(nfl,7) ! density cfs(nfl)%av(9) = cfs(nfl)%pos(1) cfs(nfl)%av(10)= cfs(nfl)%pos(2) cfs(nfl)%av(8) = cfs(nfl)%pos(3) cfs(nfl)%mirr = cmirr(nfl) cfs(nfl)%mixl = cmixl(nfl) cfs(nfl)%phy_conc = cphy_conc(nfl) cfs(nfl)%depth_max_phy_conc = cdepth_max_phy_conc(nfl) end do if (rt%iss=0.0) then alotmp=alomin alomin=alomax alomax=alotmp+360.0 elseif(flts_pos<0) then alotmp=alomax alomax=alomin alomin=alotmp-360.0 end if if (alomin>alomax) then alotmp=alomin alomin=alomax alomax=alotmp end if end subroutine correct_longitude C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAS180809 the grids is not straight north west, so we have to check if the C point is really inside the grid cell subroutine check_square(floc,inlon,inlat,is_inside) implicit none real, intent(in) :: floc(2), inlon(2,2), inlat(2,2) logical, intent(out) :: is_inside integer :: i,j,is,js,in,jn,iw,jw,ie,je real :: lon(2,2), lat(2,2) real :: tmpl real :: slope, slope1, slope2, slope3, slope4 real :: cross1, cross2, cross3, cross4 lon=inlon lat=inlat if (max(lon(1,1),lon(1,2),lon(2,1),lon(2,2)) - & min(lon(1,1),lon(1,2),lon(2,1),lon(2,2)) > 100.0) then if (floc(1) >= 0) then do i=1,2 do j=1,2 if (lon(i,j) < 0) then lon(i,j)=lon(i,j)+360; end if end do end do elseif (floc(1) < 0) then do i=1,2 do j=1,2 if (lon(i,j) >= 0) then lon(i,j)=lon(i,j)-360; end if end do end do end if end if tmpl=1000 do i=1,2 do j=1,2 if (lat(i,j)1000.0.or. & abs(min(slope1,slope2,slope3, slope4))<0.001) then is_inside=.true. return else cross1=lat(is,js)-slope1*lon(is,js); cross2=lat(is,js)-slope2*lon(is,js); cross3=lat(in,jn)-slope3*lon(in,jn); cross4=lat(in,jn)-slope4*lon(in,jn); endif if ((slope1/slope3 > 0).and.(slope2/slope4) > 0) then if (floc(2) > slope1*floc(1)+cross1) then if (floc(2) > slope2*floc(1)+cross2) then if (floc(2) < slope3*floc(1)+cross3) then if (floc(2) < slope4*floc(1)+cross4) then is_inside=.true. endif endif endif endif elseif ((slope1/slope3 < 0).and.(slope2/slope4) > 0) then if (floc(2) > slope1*floc(1)+cross1) then if (floc(2) > slope2*floc(1)+cross2) then if (floc(2) > slope3*floc(1)+cross3) then if (floc(2) < slope4*floc(1)+cross4) then is_inside=.true. endif endif endif endif elseif ((slope1/slope3 > 0).and.(slope2/slope4) < 0) then if (floc(2) > slope1*floc(1)+cross1) then if (floc(2) > slope2*floc(1)+cross2) then if (floc(2) < slope3*floc(1)+cross3) then if (floc(2) > slope4*floc(1)+cross4) then is_inside=.true. endif endif endif endif endif end subroutine check_square C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAS260606 subroutine to move floats in the vertical (mimicing vertical migration) C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine dvm(floatdepth,delt,direction,maxdepth,mindepth) real,intent(inout) :: floatdepth real,intent(in) :: delt, maxdepth, mindepth integer,intent(in) :: direction c 0=ned, 1=opp real,parameter :: speed=20.0/3600.0 if (direction==0) then floatdepth=min(floatdepth+delt*speed,maxdepth) pause else if (direction==1) then floatdepth=max(floatdepth-delt*speed,mindepth) endif end subroutine dvm C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cc c c c> Revision history: c> c> Nov 2006 - 1st module version