#include "cppdefs.h" MODULE step_floats_mod #if defined NONLINEAR && defined FLOATS ! !svn $Id: step_floats.F 1151 2023-02-09 03:08:53Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2023 The ROMS/TOMS Group John M. Klinck ! ! Licensed under a MIT/X style license ! ! See License_ROMS.txt ! !======================================================================= ! ! ! This routine time-steps simulated floats trajectories using a ! ! fourth-order Milne predictor and fourth-order Hamming corrector. ! ! ! ! Vertical diffusion is optionally represented by a random walk, ! ! in which case a forward scheme is used for vertical displacement. ! ! The probability distribution for the vertical displacement is ! ! Gaussian and includes a correction for the vertical gradient in ! ! diffusion coefficient ! ! ! # ifdef FLOAT_BIOLOGY ! If biological behavior is activated, the biological float vertical ! ! velocity is scaled to grid units by dividing by the thickness (Hz), ! ! according to the model omega vertical velocity: ! ! ! ! track(iwbio,:,:) / Hz ! ! ! ! A new field is added to the track array, track(i1ohz,:,;), having ! ! the values of 1/Hz for each float. This facilitates the scaling. ! ! ! # endif !======================================================================= ! implicit none ! PRIVATE PUBLIC :: step_floats ! CONTAINS ! !*********************************************************************** SUBROUTINE step_floats (ng, Lstr, Lend) !*********************************************************************** ! USE mod_param USE mod_floats USE mod_stepping USE mod_forces ! ! Imported variable declarations. ! integer, intent(in) :: ng, Lstr, Lend ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__ ! # ifdef PROFILE CALL wclock_on (ng, iNLM, 10, __LINE__, MyFile) # endif CALL step_floats_tile (ng, Lstr, Lend, & & knew(ng), nnew(ng), nfm3(ng), nfm2(ng), & & nfm1(ng), nf(ng), nfp1(ng), & & DRIFTER(ng) % bounded, & & DRIFTER(ng) % Ftype, & & DRIFTER(ng) % Tinfo, & & DRIFTER(ng) % Fz0, & # if defined SOLVE3D && defined FLOAT_STICKY & DRIFTER(ng) % stuck, & & DRIFTER(ng) % Ucrit, & & DRIFTER(ng) % Rho_mp, & & DRIFTER(ng) % a_axis, & & DRIFTER(ng) % b_axis, & & DRIFTER(ng) % c_axis, & # endif & DRIFTER(ng) % track) # ifdef PROFILE CALL wclock_off (ng, iNLM, 10, __LINE__, MyFile) # endif ! RETURN END SUBROUTINE step_floats ! !*********************************************************************** SUBROUTINE step_floats_tile (ng, Lstr, Lend, & & knew, nnew, & & nfm3, nfm2, nfm1, nf, nfp1, & & bounded, Ftype, Tinfo, Fz0, & # if defined SOLVE3D && defined FLOAT_STICKY & stuck, & & Ucrit, & & Rho_mp, & & a_axis, & & b_axis, & & c_axis, & # endif & track) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_floats USE mod_grid USE mod_iounits USE mod_ncparam USE mod_ocean USE mod_scalars USE mod_forces ! # ifdef FLOAT_BIOLOGY USE biology_floats_mod, ONLY : biology_floats # endif # ifdef DISTRIBUTE USE distribute_mod, ONLY : mp_collect # endif USE interp_floats_mod # if defined SOLVE3D && defined FLOAT_VWALK USE vwalk_floats_mod, ONLY : vwalk_floats # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, Lstr, Lend integer, intent(in) :: knew, nnew, nfm3, nfm2, nfm1, nf, nfp1 ! # ifdef ASSUMED_SHAPE integer, intent(in) :: Ftype(:) real(r8), intent(in) :: Tinfo(0:,:) real(r8), intent(in) :: Fz0(:) logical, intent(inout) :: bounded(:) # if defined SOLVE3D && defined FLOAT_STICKY logical, intent(inout) :: stuck(:) real(r8), intent(in) :: Ucrit(:) real(r8), intent(in) :: Rho_mp(:) real(r8), intent(in) :: a_axis(:) real(r8), intent(in) :: b_axis(:) real(r8), intent(in) :: c_axis(:) # endif real(r8), intent(inout) :: track(:,0:,:) # else integer, intent(in) :: Ftype(Nfloats(ng)) real(r8), intent(in) :: Tinfo(0:izrhs,Nfloats(ng)) real(r8), intent(in) :: Fz0(Nfloats(ng)) logical, intent(inout) :: bounded(Nfloats(ng)) # if defined SOLVE3D && defined FLOAT_STICKY logical, intent(inout) :: stuck(Nfloats(ng)) real(r8), intent(in) :: Ucrit(Nfloats(ng)) real(r8), intent(in) :: Rho_mp(Nfloats(ng)) real(r8), intent(in) :: a_axis(Nfloats(ng)) real(r8), intent(in) :: b_axis(Nfloats(ng)) real(r8), intent(in) :: c_axis(Nfloats(ng)) # endif real(r8), intent(inout) :: track(NFV(ng),0:NFT,Nfloats(ng)) # endif ! ! Local variable declarations. ! logical, parameter :: Gmask = .FALSE. # ifdef MASKING logical, parameter :: Lmask = .TRUE. # else logical, parameter :: Lmask = .FALSE. # endif logical, dimension(Lstr:Lend) :: my_thread integer :: LBi, UBi, LBj, UBj integer :: Ir, Jr, Npts, i, i1, i2, j, j1, j2, itrc, l, k, m1, m2, tmp1, tmp2 real(r8), parameter :: Fspv = 0.0_r8, nu_molecular = 1.0D-06, OneThird = 1.0_r8/3.0_r8, Two = 2.0_r8, MMtoM = 1000.0_r8 real(r8) :: cff1, cff2, cff3, cff4, cff5, cff6, cff7, cff8, cff9 real(r8) :: oHz, p1, p2, q1, q2, xrhs, yrhs, zrhs, zfloat real(r8) :: HalfDT real(r8) :: taux, tauy, bstress, temp, salt, density, ustarf, ds, rhomp, theta_cr, aaxis, baxis, caxis, Dequi real(r8), dimension(Lstr:Lend) :: nudg # ifdef DISTRIBUTE real(r8) :: Xstr, Xend, Ystr, Yend real(r8), dimension(Nfloats(ng)*NFV(ng)*(NFT+1)) :: Fwrk # endif ! ! Set tile array bounds. ! LBi=LBOUND(GRID(ng)%h,DIM=1) UBi=UBOUND(GRID(ng)%h,DIM=1) LBj=LBOUND(GRID(ng)%h,DIM=2) UBj=UBOUND(GRID(ng)%h,DIM=2) # ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! In distributed-memory configuration, determine which node bounds the ! current location of the floats. Assign unbounded floats to the master ! node. !----------------------------------------------------------------------- ! ! The strategy here is to build a switch that processes only the floats ! contained within the tile node. The trajectory data for unbounded ! floats is initialized to Fspv. These values are used during the ! collection step at the end of the routine. Since a SUM reduction is ! carried-out, setting Fspv to zero means the floats only contribute in ! their own tile. ! Npts=NFV(ng)*(NFT+1)*Nfloats(ng) Xstr=REAL(BOUNDS(ng)%Istr(MyRank),r8)-0.5_r8 Xend=REAL(BOUNDS(ng)%Iend(MyRank),r8)+0.5_r8 Ystr=REAL(BOUNDS(ng)%Jstr(MyRank),r8)-0.5_r8 Yend=REAL(BOUNDS(ng)%Jend(MyRank),r8)+0.5_r8 DO l=Lstr,Lend my_thread(l)=.FALSE. IF ((Xstr.le.track(ixgrd,nf,l)).and. & & (track(ixgrd,nf,l).lt.Xend).and. & & (Ystr.le.track(iygrd,nf,l)).and. & & (track(iygrd,nf,l).lt.Yend)) THEN my_thread(l)=.TRUE. ELSE IF (Master.and.(.not.bounded(l))) THEN my_thread(l)=.TRUE. ELSE DO j=0,NFT DO i=1,NFV(ng) track(i,j,l)=Fspv END DO END DO END IF END DO # else ! !----------------------------------------------------------------------- ! Initialize. !----------------------------------------------------------------------- ! DO l=Lstr,Lend my_thread(l)=.TRUE. END DO # endif # if !(defined SOLVE3D && defined FLOAT_VWALK) DO l=Lstr,Lend nudg(l)=0.0_r8 END DO # endif # if defined SOLVE3D && defined FLOAT_VWALK ! !----------------------------------------------------------------------- ! Compute vertical positions due to vertical random walk, predictor ! step. !----------------------------------------------------------------------- ! CALL vwalk_floats (ng, Lstr, Lend, .TRUE., my_thread, nudg) # endif ! !----------------------------------------------------------------------- ! Predictor step: compute first guess floats locations using a ! 4th-order Milne time-stepping scheme. !----------------------------------------------------------------------- ! cff1=8.0_r8/3.0_r8 cff2=4.0_r8/3.0_r8 DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN track(ixgrd,nfp1,l)=track(ixgrd,nfm3,l)+ & & dt(ng)*(cff1*track(ixrhs,nf ,l)- & & cff2*track(ixrhs,nfm1,l)+ & & cff1*track(ixrhs,nfm2,l)) track(iygrd,nfp1,l)=track(iygrd,nfm3,l)+ & & dt(ng)*(cff1*track(iyrhs,nf ,l)- & & cff2*track(iyrhs,nfm1,l)+ & & cff1*track(iyrhs,nfm2,l)) # if defined SOLVE3D && !defined FLOAT_VWALK ! ! Compute vertical position (grid units) 3D Lagrangian floats. ! IF (Ftype(l).eq.flt_Lagran) THEN # ifdef FLOAT_BIOLOGY track(izgrd,nfp1,l)=track(izgrd,nfm3,l)+ & & dt(ng)*(cff1*track(izrhs,nf ,l)- & & cff2*track(izrhs,nfm1,l)+ & & cff1*track(izrhs,nfm2,l)+ & & cff1*track(iwbio,nf ,l)* & & track(i1oHz,nf ,l)- & & cff2*track(iwbio,nfm1,l)* & & track(i1oHz,nfm1,l)+ & & cff1*track(iwbio,nfm2,l)* & & track(i1oHz,nfm2,l)) # else track(izgrd,nfp1,l)=track(izgrd,nfm3,l)+ & & dt(ng)*(cff1*track(izrhs,nf ,l)- & & cff2*track(izrhs,nfm1,l)+ & & cff1*track(izrhs,nfm2,l)) # endif ! ! Compute vertical position (grid units) for isobaric floats ! (p=g*(z+zeta)=constant) or geopotential floats (constant depth). ! Use bilinear interpolation to determine vertical position. ! ELSE IF ((Ftype(l).eq.flt_Isobar).or. & & (Ftype(l).eq.flt_Geopot)) THEN Ir=INT(track(ixgrd,nfp1,l)) Jr=INT(track(iygrd,nfp1,l)) ! i1=MIN(MAX(Ir ,0),Lm(ng)+1) i2=MIN(MAX(Ir+1,1),Lm(ng)+1) j1=MIN(MAX(Jr ,0),Mm(ng)+1) j2=MIN(MAX(Jr+1,0),Mm(ng)+1) ! p2=REAL(i2-i1,r8)*(track(ixgrd,nfp1,l)-REAL(i1,r8)) q2=REAL(j2-j1,r8)*(track(iygrd,nfp1,l)-REAL(j1,r8)) p1=1.0_r8-p2 q1=1.0_r8-q2 # ifdef MASKING cff7=p1*q1*GRID(ng)%z_w(i1,j1,N(ng))*GRID(ng)%rmask(i1,j1)+ & & p2*q1*GRID(ng)%z_w(i2,j1,N(ng))*GRID(ng)%rmask(i2,j1)+ & & p1*q2*GRID(ng)%z_w(i1,j2,N(ng))*GRID(ng)%rmask(i1,j2)+ & & p2*q2*GRID(ng)%z_w(i2,j2,N(ng))*GRID(ng)%rmask(i2,j2) cff8=p1*q1*GRID(ng)%rmask(i1,j1)+ & & p2*q1*GRID(ng)%rmask(i2,j1)+ & & p1*q2*GRID(ng)%rmask(i1,j2)+ & & p2*q2*GRID(ng)%rmask(i2,j2) cff9=0.0_r8 IF (cff8.gt.0.0_r8) cff9=cff7/cff8 # else cff9=p1*q1*GRID(ng)%z_w(i1,j1,N(ng))+ & & p2*q1*GRID(ng)%z_w(i2,j1,N(ng))+ & & p1*q2*GRID(ng)%z_w(i1,j2,N(ng))+ & & p2*q2*GRID(ng)%z_w(i2,j2,N(ng)) # endif cff6=cff9 ! IF (Ftype(l).eq.flt_Geopot) THEN zfloat=Fz0(l) ELSE IF (Ftype(l).eq.flt_Isobar) THEN zfloat=Fz0(l)+cff9 END IF ! DO k=N(ng)-1,0,-1 # ifdef MASKING cff7=p1*q1*GRID(ng)%z_w(i1,j1,k)*GRID(ng)%rmask(i1,j1)+ & & p2*q1*GRID(ng)%z_w(i2,j1,k)*GRID(ng)%rmask(i2,j1)+ & & p1*q2*GRID(ng)%z_w(i1,j2,k)*GRID(ng)%rmask(i1,j2)+ & & p2*q2*GRID(ng)%z_w(i2,j2,k)*GRID(ng)%rmask(i2,j2) cff8=p1*q1*GRID(ng)%rmask(i1,j1)+ & & p2*q1*GRID(ng)%rmask(i2,j1)+ & & p1*q2*GRID(ng)%rmask(i1,j2)+ & & p2*q2*GRID(ng)%rmask(i2,j2) IF (cff8.gt.0.0_r8) THEN cff5=cff7/cff8 ELSE cff5=0.0_r8 END IF # else cff5=p1*q1*GRID(ng)%z_w(i1,j1,k)+ & & p2*q1*GRID(ng)%z_w(i2,j1,k)+ & & p1*q2*GRID(ng)%z_w(i1,j2,k)+ & & p2*q2*GRID(ng)%z_w(i2,j2,k) # endif IF ((zfloat-cff5)*(cff6-zfloat).ge.0.0_r8) THEN track(izgrd,nfp1,l)=REAL(k,r8)+(zfloat-cff5)/(cff6-cff5) END IF cff6=cff5 END DO END IF # endif END IF END DO ! !----------------------------------------------------------------------- ! Calculate slopes at new time-step. !----------------------------------------------------------------------- ! # ifdef SOLVE3D CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, ixrhs, isUvel, & & -u3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % u(:,:,:,nnew), & & my_thread, bounded, track) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, iyrhs, isVvel, & & -v3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % v(:,:,:,nnew), & & my_thread, bounded, track) # if !defined FLOAT_VWALK CALL interp_floats (ng, LBi, UBi, LBj, UBj, 0, N(ng), & & Lstr, Lend, nfp1, izrhs, isBw3d, & & -w3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % W, & & my_thread, bounded, track) # endif # if defined FLOAT_BIOLOGY CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, i1oHz, isBr3d, & & r3dvar, Gmask, Fspv, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & GRID(ng) % Hz, & & my_thread, bounded, track) DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN track(i1oHz,nfp1,l)=1.0_r8/track(i1oHz,nfp1,l) END IF END DO # endif # else CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, 1, & & Lstr, Lend, nfp1, ixrhs, isUbar, & & -u2dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % ubar(:,:,knew), & & my_thread, bounded, track) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, 1, & & Lstr, Lend, nfp1, iyrhs, isVbar, & & -v2dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % vbar(:,:,knew), & & my_thread, bounded, track) # endif # ifdef FLOAT_BIOLOGY ! !----------------------------------------------------------------------- ! Compute biological float behavior, predictor step. !----------------------------------------------------------------------- ! ! Interpolate tracer to the predictor step locations. These values ! are used in the "biology_floats" routine. ! DO itrc=1,NT(ng) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, ifTvar(itrc), & & isTvar(itrc), r3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % t(:,:,:,nnew,itrc), & & my_thread, bounded, track) END DO ! ! Biological behavior predictor step. ! CALL biology_floats (ng, Lstr, Lend, .TRUE., my_thread) # endif ! !----------------------------------------------------------------------- ! Corrector step: correct floats locations using a 4th-order ! Hamming time-stepping scheme. !----------------------------------------------------------------------- ! cff1=9.0_r8/8.0_r8 cff2=1.0_r8/8.0_r8 cff3=3.0_r8/8.0_r8 cff4=6.0_r8/8.0_r8 DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN track(ixgrd,nfp1,l)=cff1*track(ixgrd,nf ,l)- & & cff2*track(ixgrd,nfm2,l)+ & & dt(ng)*(cff3*track(ixrhs,nfp1,l)+ & & cff4*track(ixrhs,nf ,l)- & & cff3*track(ixrhs,nfm1,l)) track(iygrd,nfp1,l)=cff1*track(iygrd,nf ,l)- & & cff2*track(iygrd,nfm2,l)+ & & dt(ng)*(cff3*track(iyrhs,nfp1,l)+ & & cff4*track(iyrhs,nf ,l)- & & cff3*track(iyrhs,nfm1,l)) # if defined SOLVE3D && !defined FLOAT_VWALK ! ! Compute vertical position (grid units) 3D Lagrangian floats. ! IF (Ftype(l).eq.flt_Lagran) THEN # if defined FLOAT_BIOLOGY track(izgrd,nfp1,l)=cff1*track(izgrd,nf ,l)- & & cff2*track(izgrd,nfm2,l)+ & & dt(ng)*(cff3*track(izrhs,nfp1,l)+ & & cff4*track(izrhs,nf ,l)- & & cff3*track(izrhs,nfm1,l)+ & & cff3*track(iwbio,nfp1,l)* & & track(i1oHz,nfp1,l)+ & & cff4*track(iwbio,nf ,l)* & & track(i1oHz,nf ,l)- & & cff3*track(iwbio,nfm1,l)* & & track(i1oHz,nf ,l)) # else track(izgrd,nfp1,l)=cff1*track(izgrd,nf ,l)- & & cff2*track(izgrd,nfm2,l)+ & & dt(ng)*(cff3*track(izrhs,nfp1,l)+ & & cff4*track(izrhs,nf ,l)- & & cff3*track(izrhs,nfm1,l)) # endif ! ! Compute vertical position (grid units) for isobaric floats ! (p=g*(z+zeta)=constant) or geopotential floats (constant depth). ! Use bilinear interpolation to determine vertical position. ! ELSE IF ((Ftype(l).eq.flt_Isobar).or. & & (Ftype(l).eq.flt_Geopot)) THEN Ir=INT(track(ixgrd,nfp1,l)) Jr=INT(track(iygrd,nfp1,l)) ! i1=MIN(MAX(Ir ,0),Lm(ng)+1) i2=MIN(MAX(Ir+1,1),Lm(ng)+1) j1=MIN(MAX(Jr ,0),Mm(ng)+1) j2=MIN(MAX(Jr+1,0),Mm(ng)+1) ! p2=REAL(i2-i1,r8)*(track(ixgrd,nfp1,l)-REAL(i1,r8)) q2=REAL(j2-j1,r8)*(track(iygrd,nfp1,l)-REAL(j1,r8)) p1=1.0_r8-p2 q1=1.0_r8-q2 # ifdef MASKING cff7=p1*q1*GRID(ng)%z_w(i1,j1,N(ng))*GRID(ng)%rmask(i1,j1)+ & & p2*q1*GRID(ng)%z_w(i2,j1,N(ng))*GRID(ng)%rmask(i2,j1)+ & & p1*q2*GRID(ng)%z_w(i1,j2,N(ng))*GRID(ng)%rmask(i1,j2)+ & & p2*q2*GRID(ng)%z_w(i2,j2,N(ng))*GRID(ng)%rmask(i2,j2) cff8=p1*q1*GRID(ng)%rmask(i1,j1)+ & & p2*q1*GRID(ng)%rmask(i2,j1)+ & & p1*q2*GRID(ng)%rmask(i1,j2)+ & & p2*q2*GRID(ng)%rmask(i2,j2) IF (cff8.gt.0.0_r8) THEN cff9=cff7/cff8 ELSE cff9=0.0_r8 END IF # else cff9=p1*q1*GRID(ng)%z_w(i1,j1,N(ng))+ & & p2*q1*GRID(ng)%z_w(i2,j1,N(ng))+ & & p1*q2*GRID(ng)%z_w(i1,j2,N(ng))+ & & p2*q2*GRID(ng)%z_w(i2,j2,N(ng)) # endif cff6=cff9 ! IF (Ftype(l).eq.flt_Geopot) THEN zfloat=Fz0(l) ELSE IF (Ftype(l).eq.flt_Isobar) THEN zfloat=Fz0(l)+cff9 END IF ! DO k=N(ng)-1,0,-1 # ifdef MASKING cff7=p1*q1*GRID(ng)%z_w(i1,j1,k)*GRID(ng)%rmask(i1,j1)+ & & p2*q1*GRID(ng)%z_w(i2,j1,k)*GRID(ng)%rmask(i2,j1)+ & & p1*q2*GRID(ng)%z_w(i1,j2,k)*GRID(ng)%rmask(i1,j2)+ & & p2*q2*GRID(ng)%z_w(i2,j2,k)*GRID(ng)%rmask(i2,j2) cff8=p1*q1*GRID(ng)%rmask(i1,j1)+ & & p2*q1*GRID(ng)%rmask(i2,j1)+ & & p1*q2*GRID(ng)%rmask(i1,j2)+ & & p2*q2*GRID(ng)%rmask(i2,j2) cff5=0.0_r8 IF (cff8.gt.0.0_r8) cff5=cff7/cff8 # else cff5=p1*q1*GRID(ng)%z_w(i1,j1,k)+ & & p2*q1*GRID(ng)%z_w(i2,j1,k)+ & & p1*q2*GRID(ng)%z_w(i1,j2,k)+ & & p2*q2*GRID(ng)%z_w(i2,j2,k) # endif IF ((zfloat-cff5)*(cff6-zfloat).ge.0.0_r8) THEN track(izgrd,nfp1,l)=REAL(k,r8)+(zfloat-cff5)/(cff6-cff5) END IF cff6=cff5 END DO END IF # endif END IF END DO ! !----------------------------------------------------------------------- ! Determine floats status. !----------------------------------------------------------------------- ! IF (EWperiodic(ng)) THEN cff1=REAL(Lm(ng),r8) DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN IF (track(ixgrd,nfp1,l).ge.REAL(Lm(ng)+1,r8)-0.5_r8) THEN track(ixgrd,nfp1,l)=track(ixgrd,nfp1,l)-cff1 track(ixgrd,nf ,l)=track(ixgrd,nf ,l)-cff1 track(ixgrd,nfm1,l)=track(ixgrd,nfm1,l)-cff1 track(ixgrd,nfm2,l)=track(ixgrd,nfm2,l)-cff1 track(ixgrd,nfm3,l)=track(ixgrd,nfm3,l)-cff1 ELSE IF (track(ixgrd,nfp1,l).lt.0.5_r8) THEN track(ixgrd,nfp1,l)=cff1+track(ixgrd,nfp1,l) track(ixgrd,nf ,l)=cff1+track(ixgrd,nf ,l) track(ixgrd,nfm1,l)=cff1+track(ixgrd,nfm1,l) track(ixgrd,nfm2,l)=cff1+track(ixgrd,nfm2,l) track(ixgrd,nfm3,l)=cff1+track(ixgrd,nfm3,l) END IF END IF END DO # ifdef DISTRIBUTE IF (NtileI(ng).gt.1) THEN Fwrk=RESHAPE(track,(/Npts/)) CALL mp_collect (ng, iNLM, Npts, Fspv, Fwrk) track=RESHAPE(Fwrk,(/NFV(ng),NFT+1,Nfloats(ng)/)) DO l=Lstr,Lend IF ((Xstr.le.track(ixgrd,nfp1,l)).and. & & (track(ixgrd,nfp1,l).lt.Xend).and. & & (Ystr.le.track(iygrd,nfp1,l)).and. & & (track(iygrd,nfp1,l).lt.Yend)) THEN my_thread(l)=.TRUE. ELSE IF (Master.and.(.not.bounded(l))) THEN my_thread(l)=.TRUE. ELSE my_thread(l)=.FALSE. DO j=0,NFT DO i=1,NFV(ng) track(i,j,l)=Fspv END DO END DO END IF END DO END IF # endif ELSE DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN IF ((track(ixgrd,nfp1,l).ge.REAL(Lm(ng)+1,r8)-0.5_r8).or. & & (track(ixgrd,nfp1,l).lt.0.5_r8)) THEN bounded(l)=.FALSE. END IF END IF END DO END IF ! IF (NSperiodic(ng)) THEN cff1=REAL(Mm(ng),r8) DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN IF (track(iygrd,nfp1,l).ge.REAL(Mm(ng)+1,r8)-0.5_r8) THEN track(iygrd,nfp1,l)=track(iygrd,nfp1,l)-cff1 track(iygrd,nf ,l)=track(iygrd,nf ,l)-cff1 track(iygrd,nfm1,l)=track(iygrd,nfm1,l)-cff1 track(iygrd,nfm2,l)=track(iygrd,nfm2,l)-cff1 track(iygrd,nfm3,l)=track(iygrd,nfm3,l)-cff1 ELSE IF (track(iygrd,nfp1,l).lt.0.5_r8) THEN track(iygrd,nfp1,l)=cff1+track(iygrd,nfp1,l) track(iygrd,nf ,l)=cff1+track(iygrd,nf ,l) track(iygrd,nfm1,l)=cff1+track(iygrd,nfm1,l) track(iygrd,nfm2,l)=cff1+track(iygrd,nfm2,l) track(iygrd,nfm3,l)=cff1+track(iygrd,nfm3,l) END IF END IF END DO # ifdef DISTRIBUTE IF (NtileJ(ng).gt.1) THEN Fwrk=RESHAPE(track,(/Npts/)) CALL mp_collect (ng, iNLM, Npts, Fspv, Fwrk) track=RESHAPE(Fwrk,(/NFV(ng),NFT+1,Nfloats(ng)/)) DO l=Lstr,Lend IF ((Xstr.le.track(ixgrd,nfp1,l)).and. & & (track(ixgrd,nfp1,l).lt.Xend).and. & & (Ystr.le.track(iygrd,nfp1,l)).and. & & (track(iygrd,nfp1,l).lt.Yend)) THEN my_thread(l)=.TRUE. ELSE IF (Master.and.(.not.bounded(l))) THEN my_thread(l)=.TRUE. ELSE my_thread(l)=.FALSE. DO j=0,NFT DO i=1,NFV(ng) track(i,j,l)=Fspv END DO END DO END IF END DO END IF # endif ELSE DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN IF ((track(iygrd,nfp1,l).ge.REAL(Mm(ng)+1,r8)-0.5_r8).or. & & (track(iygrd,nfp1,l).lt.0.5_r8)) THEN bounded(l)=.FALSE. END IF END IF END DO END IF ! !----------------------------------------------------------------------- ! If appropriate, activate the release of new floats and set initial ! positions for all time levels. !----------------------------------------------------------------------- ! HalfDT=0.5_r8*dt(ng) DO l=Lstr,Lend IF (.not.bounded(l).and. & & (time(ng)-HalfDT.le.Tinfo(itstr,l).and. & & time(ng)+HalfDT.gt.Tinfo(itstr,l))) THEN bounded(l)=.TRUE. IF ((Tinfo(ixgrd,l).lt.0.5_r8).or. & & (Tinfo(iygrd,l).lt.0.5_r8).or. & & (Tinfo(ixgrd,l).gt.REAL(Lm(ng),r8)+0.5_r8).or. & & (Tinfo(iygrd,l).gt.REAL(Mm(ng),r8)+0.5_r8)) THEN bounded(l)=.FALSE. ! outside application grid END IF # if defined SOLVE3D && defined FLOAT_STICKY stuck(l)=.FALSE. # endif # ifdef DISTRIBUTE IF ((Xstr.le.Tinfo(ixgrd,l)).and. & & (Tinfo(ixgrd,l).lt.Xend).and. & & (Ystr.le.Tinfo(iygrd,l)).and. & & (Tinfo(iygrd,l).lt.Yend).and.bounded(l)) THEN DO j=0,NFT track(ixgrd,j,l)=Tinfo(ixgrd,l) track(iygrd,j,l)=Tinfo(iygrd,l) # ifdef SOLVE3D track(izgrd,j,l)=Tinfo(izgrd,l) # endif END DO my_thread(l)=.TRUE. ELSE my_thread(l)=.FALSE. DO j=0,NFT DO i=1,NFV(ng) track(i,j,l)=Fspv END DO END DO END IF # else IF (bounded(l)) THEN my_thread(l)=.TRUE. DO j=0,NFT track(ixgrd,j,l)=Tinfo(ixgrd,l) track(iygrd,j,l)=Tinfo(iygrd,l) # ifdef SOLVE3D track(izgrd,j,l)=Tinfo(izgrd,l) # endif END DO END IF # endif END IF END DO ! !----------------------------------------------------------------------- ! Calculate slopes with corrected locations. !----------------------------------------------------------------------- ! # ifdef SOLVE3D CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, ixrhs, isUvel, & & -u3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % u(:,:,:,nnew), & & my_thread, bounded, track) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, iyrhs, isVvel, & & -v3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % v(:,:,:,nnew), & & my_thread, bounded, track) # if !defined FLOAT_VWALK CALL interp_floats (ng, LBi, UBi, LBj, UBj, 0, N(ng), & & Lstr, Lend, nfp1, izrhs, isBw3d, & & -w3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % W, & & my_thread, bounded, track) # endif # ifdef FLOAT_BIOLOGY CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, i1oHz, isBr3d, & & r3dvar, Gmask, Fspv, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & GRID(ng) % Hz, & & my_thread, bounded, track) DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN track(i1oHz,nfp1,l)=1.0_r8/track(i1oHz,nfp1,l) END IF END DO # endif # else CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, 1, & & Lstr, Lend, nfp1, ixrhs, isUbar, & & -u2dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % ubar(:,:,knew), & & my_thread, bounded, track) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, 1, & & Lstr, Lend, nfp1, iyrhs, isVbar, & & -v2dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % vbar(:,:,knew), & & my_thread, bounded, track) # endif ! ! If newly released floats, initialize slopes at all time levels. ! DO l=Lstr,Lend IF (my_thread(l).and.bounded(l).and. & & (time(ng)-HalfDT.le.Tinfo(itstr,l).and. & & time(ng)+HalfDT.gt.Tinfo(itstr,l))) THEN xrhs=track(ixrhs,nfp1,l) yrhs=track(iyrhs,nfp1,l) # ifdef SOLVE3D zrhs=track(izrhs,nfp1,l) # endif # ifdef FLOAT_BIOLOGY oHz=track(i1oHz,nfp1,l) # endif DO i=0,NFT track(ixrhs,i,l)=xrhs track(iyrhs,i,l)=yrhs # ifdef SOLVE3D track(izrhs,i,l)=zrhs # endif # ifdef FLOAT_BIOLOGY track(i1oHz,i,l)=oHz # endif END DO END IF END DO ! !----------------------------------------------------------------------- ! Interpolate various output variables at the corrected locations. !----------------------------------------------------------------------- ! IF (spherical) THEN CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, 1, & & Lstr, Lend, nfp1, iflon, isBr2d, & & r2dvar, Gmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef SOLVE3D & GRID(ng) % Hz, & # endif # ifdef MASKING & GRID(ng) % rmask, & # endif & GRID(ng) % lonr, & & my_thread, bounded, track) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, 1, & & Lstr, Lend, nfp1, iflat, isBr2d, & & r2dvar, Gmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef SOLVE3D & GRID(ng) % Hz, & # endif # ifdef MASKING & GRID(ng) % rmask, & # endif & GRID(ng) % latr, & & my_thread, bounded, track) ELSE CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, 1, & & Lstr, Lend, nfp1, iflon, isBr2d, & & r2dvar, Gmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef SOLVE3D & GRID(ng) % Hz, & # endif # ifdef MASKING & GRID(ng) % rmask, & # endif & GRID(ng) % xr, & & my_thread, bounded, track) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, 1, & & Lstr, Lend, nfp1, iflat, isBr2d, & & r2dvar, Gmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef SOLVE3D & GRID(ng) % Hz, & # endif # ifdef MASKING & GRID(ng) % rmask, & # endif & GRID(ng) % yr, & & my_thread, bounded, track) END IF # ifdef SOLVE3D CALL interp_floats (ng, LBi, UBi, LBj, UBj, 0, N(ng), & & Lstr, Lend, nfp1, idpth, isBw3d, & & w3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & GRID(ng) % z_w, & & my_thread, bounded, track) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, ifden, isBr3d, & & r3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % rho, & & my_thread, bounded, track) DO itrc=1,NT(ng) CALL interp_floats (ng, LBi, UBi, LBj, UBj, 1, N(ng), & & Lstr, Lend, nfp1, ifTvar(itrc), & & isTvar(itrc), r3dvar, Lmask, spval, nudg, & & GRID(ng) % pm, & & GRID(ng) % pn, & & GRID(ng) % Hz, & # ifdef MASKING & GRID(ng) % rmask, & # endif & OCEAN(ng) % t(:,:,:,nnew,itrc), & & my_thread, bounded, track) END DO # endif # ifdef FLOAT_BIOLOGY ! !----------------------------------------------------------------------- ! Compute biological float behavior, corrector step. !----------------------------------------------------------------------- ! CALL biology_floats (ng, Lstr, Lend, .FALSE., my_thread) # endif # if defined SOLVE3D && defined FLOAT_VWALK && !defined VWALK_FORWARD ! !----------------------------------------------------------------------- ! Compute vertical positions due to vertical random walk, corrector ! step. !----------------------------------------------------------------------- ! CALL vwalk_floats (ng, Lstr, Lend, .FALSE., my_thread, nudg) # endif # ifdef SOLVE3D # ifdef FLOAT_STICKY ! ! Floats that hit the surface are reflected; floats that hit the bottom ! get stuck. ! DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN ! ! IF (track(izgrd,nfp1,l).gt.REAL(N(ng),r8)) THEN ! DO j=0,NFT ! track(izgrd,j,l)=2.0_r8*REAL(N(ng),r8)- & ! & track(izgrd,j,l) ! END DO ! ELSE IF (track(izgrd,nfp1,l).lt.0.0_r8) THEN ! DO j=0,NFT ! track(izgrd,j,l)=0.0_r8 ! END DO ! track(izgrd,nfp1,l) = 0.0_r8 ! stuck(l)=.TRUE. ! END IF ! IF (stuck(l)) THEN track(ixgrd,nfp1,l)=track(ixgrd,nf,l) track(iygrd,nfp1,l)=track(iygrd,nf,l) track(izgrd,nfp1,l)=track(izgrd,nf,l) ELSE IF (track(izgrd,nfp1,l).gt.REAL(N(ng),r8)) THEN DO j=0,NFT track(izgrd,j,l)=2.0_r8*REAL(N(ng),r8)- & & track(izgrd,j,l) END DO ELSE IF (track(izgrd,nfp1,l).lt.0.0_r8) THEN DO j=0,NFT track(izgrd,j,l)=0.0_r8 END DO ! track(izgrd,nfp1,l) = 0.0_r8 stuck(l)=.TRUE. END IF END IF ! IF (stuck(l)) ! IF (stuck(l)) THEN ! Re-label floats (stuck) according to sqrt(tau_b/rho) > or < Ucrit; for sinking ones, they are ALWAYS modelled as "oysters"... ! Step 1: compute density at the given location temp=track(ifTvar(itemp),nfp1,l) salt=track(ifTvar(isalt),nfp1,l) density = rho0 * (1.0_r8 - Tcoef(ng)*(temp-T0(ng)) + Scoef(ng)*(salt-S0(ng))) ! Step 2: compute theta_cr on the seabed rhomp = Rho_mp(l) aaxis = a_axis(l)/MMtoM baxis = b_axis(l)/MMtoM caxis = c_axis(l)/MMtoM Dequi = Two*(aaxis*baxis*caxis)**(OneThird) ds = (ABS(rhomp/density-1)*g*nu_molecular**(-Two))**(OneThird)*Dequi IF ((ds .GT. 10.0_r8) .and. (ds .LE. 20.0_r8)) THEN theta_cr = 0.04_r8*ds**(-0.1_r8) ELSE IF ((ds .GT. 4.0_r8) .and. (ds .LE. 10.0_r8)) THEN theta_cr = 0.14_r8*ds**(-0.64_r8) ELSE IF (ds .LE. 4.0_r8) THEN theta_cr = 0.24_r8*ds**(-1.0_r8) ELSE IF (ds .GT. 20.0_r8) THEN theta_cr = 0.013_r8*ds**(0.29_r8); END IF IF (Ucrit(l) .GE. 999.0_r8) THEN theta_cr = 999.0_r8 END IF theta_cr = 0.0_r8 ! Step 3: compute bottom stress i1 = MIN(INT(track(ixgrd,nfp1,1)+0.5_r8),Lm(ng)+1) i2 = MAX(i1 + 1,0) !MIN(MAX(0,INT(track(ixgrd,nfp1,l))),Lm(ng)+1) j1 = MIN(INT(track(iygrd,nfp1,1)+0.5_r8),Mm(ng)+1) j2 = MAX(j1 + 1,0) taux=((i2-0.5_r8)-track(ixgrd,nfp1,l))* & & REAL(GRID(ng)%umask(i1,j1),r8)* & & FORCES(ng)%bustr(i1,j1)+ & & (track(ixgrd,nfp1,l)-(i1-0.5_r8))* & & REAL(GRID(ng)%umask(i2,j1),r8)* & & FORCES(ng)%bustr(i2,j1) ! tauy=((j2-0.5_r8)-track(iygrd,nfp1,l))* & & REAL(GRID(ng)%vmask(i1,j1),r8)* & & FORCES(ng)%bvstr(i1,j1)+ & & (track(iygrd,nfp1,l)-(j1-0.5_r8))* & & REAL(GRID(ng)%vmask(i1,j2),r8)* & & FORCES(ng)%bvstr(i1,j2) bstress = rho0*SQRT(taux**Two + tauy**Two) ! "ustarf" is NOT friction velocity; it is bstress/((rhomp-rhow)*g*Dequi) ustarf = bstress/((rhomp-density)*g*Dequi) ! ! STEP 4: Re-labelling: IF ( (ustarf .gt. theta_cr) .or. (rhomp .lt. density) ) THEN stuck(l) = .FALSE. ! STEP 5: Re-intialise the rhs of the re-mobilised particle with the most up-to-date ! hydrodynamic forcing, as if it were newly released from this location xrhs=track(ixrhs,nfp1,l) yrhs=track(iyrhs,nfp1,l) # ifdef SOLVE3D zrhs=track(izrhs,nfp1,l) # endif # ifdef FLOAT_BIOLOGY oHz=track(i1oHz,nfp1,l) # endif DO i=0,NFT track(ixrhs,i,l)=xrhs track(iyrhs,i,l)=yrhs # ifdef SOLVE3D track(izrhs,i,l)=zrhs # endif # ifdef FLOAT_BIOLOGY track(i1oHz,i,l)=oHz # endif END DO END IF ! if ustarf > ucrit END IF !IF (stuck(l)) ! END IF END DO # else ! ! Floats that hit the surface or bottom are reflected ! DO l=Lstr,Lend IF (my_thread(l).and.bounded(l)) THEN IF (track(izgrd,nfp1,l).gt.REAL(N(ng),r8)) THEN DO j=0,NFT track(izgrd,j,l)=2.0_r8*REAL(N(ng),r8)-track(izgrd,j,l) END DO ELSE IF (track(izgrd,nfp1,l).lt.0.0_r8) THEN DO j=0,NFT track(izgrd,:,l)=-track(izgrd,:,l) END DO END IF END IF END DO # endif # endif # ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! Collect floats on all nodes. !----------------------------------------------------------------------- ! Fwrk=RESHAPE(track,(/Npts/)) CALL mp_collect (ng, iNLM, Npts, Fspv, Fwrk) track=RESHAPE(Fwrk,(/NFV(ng),NFT+1,Nfloats(ng)/)) ! ! Collect the bounded status switch. ! Fwrk=Fspv DO l=1,Nfloats(ng) IF (bounded(l)) THEN Fwrk(l)=1.0_r8 END IF END DO CALL mp_collect (ng, iNLM, Nfloats(ng), Fspv, Fwrk) DO l=1,Nfloats(ng) IF (Fwrk(l).ne.Fspv) THEN bounded(l)=.TRUE. ELSE bounded(l)=.FALSE. END IF END DO # endif ! RETURN END SUBROUTINE step_floats_tile #endif END MODULE step_floats_mod