I redid the test plots for comparison. Also, I added a check to diffu…

…sion.f95 that checks to see if the drifter is on land, and if its new displacement puts it on land, then it calculates a new displacement. With this change, the different subgrid turbulence parameterizations have consistent behavior at the coastline, which is to disallow drifters beaching.

cross.f95

@@ -72,6 +72,9 @@ subroutine cross(ijk,ia,ja,ka,r0,sp,sn,rr,uflux,vflux,wflux,ff,km,jmt,imt,do3d,d
     uu=(rg*uflux(ia,ja,ka,nsp)+rr*uflux(ia,ja,ka,nsm))*ff ! this is interpolation between time fields
     um=(rg*uflux(im,ja,ka,nsp)+rr*uflux(im,ja,ka,nsm))*ff
 
+!     print *,'x: before adding in turb vels'
+!     print *,'uu=',uu,' um=',um
+
     if(doturb==1) then   
         if(r0.ne.dble(ii)) then
             uu=uu+upr(1,2)  
@@ -85,11 +88,17 @@ subroutine cross(ijk,ia,ja,ka,r0,sp,sn,rr,uflux,vflux,wflux,ff,km,jmt,imt,do3d,d
         endif
     endif
 
+!     print *,'x: after adding in turb vels'
+!     print *,'uu=',uu,' um=',um
+
 else if(ijk.eq.2) then
     ii=ja
     uu=(rg*vflux(ia,ja  ,ka,nsp)+rr*vflux(ia,ja  ,ka,nsm))*ff
     um=(rg*vflux(ia,ja-1,ka,nsp)+rr*vflux(ia,ja-1,ka,nsm))*ff
 
+!     print *,'y: before adding in turb vels'
+!     print *,'uu=',uu,' um=',um
+
     if(doturb==1) then   
         if(r0.ne.dble(ja  )) then
             uu=uu+upr(3,2)  
@@ -103,6 +112,9 @@ subroutine cross(ijk,ia,ja,ka,r0,sp,sn,rr,uflux,vflux,wflux,ff,km,jmt,imt,do3d,d
         endif
     endif
 
+!     print *,'y: after adding in turb vels'
+!     print *,'uu=',uu,' um=',um
+
 else if(ijk.eq.3) then
     ii=ka
 ! #ifdef full_wflux

diffusion.f95

@@ -46,16 +46,21 @@ subroutine diffuse(x1, y1, z1, ib, jb, kb, dt,imt,jmt,km,kmt,dxv,dyu,dzt,h,ah,av
 real(kind=8)                                                :: xd, yd, zd, tmpX, tmpY, tmpZ
 logical                                                     :: tryAgain       
 
+! print *,''
+! print *,'start of diffusion'
+! print *,'x1=',x1,' y1=',y1,' ib=',ib,' jb=',jb
+
 tryAgain = .FALSE.
 
 ! Is particle within model area?
 ! KMT: I do not understand how the vertical checks are correct here, so I am changing them
 
 !   if(ib>=1 .AND. ib<=imt .AND. jb>=1 .AND. jb<=jmt .AND. km+1-kmt(ib,jb)<=kb .AND. kb>=1 ) then
-if(ib>=1 .AND. ib<=imt .AND. jb>=1 .AND. jb<=jmt .AND. km>=kb .AND. kb>=1 ) then
+if(ib>1 .AND. ib<imt .AND. jb>1 .AND. jb<jmt .AND. km>=kb .AND. kb>=1 ) then
     tryAgain = .TRUE.
 else
-    print *,'outside model domain in diffusion',ib,jb,km+1-kmt(ib,jb),kb
+    print *,'outside model domain in diffusion',ib,jb,km,kb
+!     print *,'outside model domain in diffusion',ib,jb,km+1-kmt(ib,jb),kb
 end if
 
 itno=0
@@ -71,6 +76,8 @@ subroutine diffuse(x1, y1, z1, ib, jb, kb, dt,imt,jmt,km,kmt,dxv,dyu,dzt,h,ah,av
     xd = xd/dxv(ib,jb)
     yd = yd/dyu(ib,jb)
     zd = zd/dzt(ib,jb,kb,1) !KMT: this should be better than dz()
+!     print *,'dzt(ib,jb,kb,1)=',dzt(ib,jb,kb,1)
+
     ! TO DO: improve so that dzt is calculated including free surface and for correct time
 !     zd = zd/dz(kb) !should be replaced for bottom box and include ssh (original note)
 
@@ -84,21 +91,31 @@ subroutine diffuse(x1, y1, z1, ib, jb, kb, dt,imt,jmt,km,kmt,dxv,dyu,dzt,h,ah,av
     tmpi = int(tmpX) + 1
     tmpj = int(tmpY) + 1
     tmpk = int(tmpZ) + 1
+!     print *,'temp dzt(ib,jb,kb,1)=',dzt(tmpi,tmpj,tmpk,1)
 
-    ! Check if particle is on an open boundary. Using rco example.
-    if(tmpi==1 .AND. tmpj>=1 .AND. tmpj<=jmt .AND. km+1-kmt(tmpi,tmpj)<=tmpk .AND. tmpk>=1 ) then
-        tryAgain = .FALSE.
-    end if
+!     ! Check if particle is on an open boundary. Using rco example.
+!     if(tmpi==1 .AND. tmpj>=1 .AND. tmpj<=jmt .AND. km+1-kmt(tmpi,tmpj)<=tmpk .AND. tmpk>=1 ) then
+!         tryAgain = .FALSE.
+!     end if
 
     ! check that column is deep enough  
     ! KMT: this is the same check as above that looks wrong and I changed
 !     if( 1<=tmpi .AND. tmpi<=imt .AND. 1<=tmpj .AND. tmpj<=jmt .AND. km+1-kmt(tmpi,tmpj)<=tmpk .AND. tmpk>=1 ) then
-    if(tmpi>=1 .AND. tmpi<=imt .AND. tmpj>=1 .AND. tmpj<=jmt .AND. km>=tmpk .AND. tmpk>=1 ) then
+    if(tmpi>1 .AND. tmpi<imt .AND. tmpj>1 .AND. tmpj<jmt .AND. km>=tmpk .AND. tmpk>=1 ) then
 !         print *,'km=',km,' km+1-kmt(tmpi,tmpj)=',km+1-kmt(tmpi,tmpj),' tmpk=',tmpk
         tryAgain = .FALSE. 
+!         print *,'found new position for drifter that is within model domain, so move on.'
         ! if false then a new position for the particle has been found and we exit the loop
     end if 
-
+
+    ! KMT: Need to check to see if new position is on masked land.
+    ! Can do this by checking dzt at the new position, which is zero when on land.
+    ! If we are on land, keep tryAgain = .TRUE. so that this displacement choice is
+    ! not used. This can alter the value of tryAgain just given above.
+    if(dzt(tmpi,tmpj,tmpk,1)==0.) then
+        tryAgain = .TRUE.
+    endif
+
     ! If tryAgain is still true, the particle is outside model area. The
     ! displacement is not saved, but we make a new try to displace.
 
@@ -120,6 +137,9 @@ subroutine diffuse(x1, y1, z1, ib, jb, kb, dt,imt,jmt,km,kmt,dxv,dyu,dzt,h,ah,av
 ib = tmpi
 jb = tmpj
 
+! print *,'After diffusion'
+! print *,'x1=',x1,' y1=',y1,' ib=',ib,' jb=',jb
+
 if(do3d==1) then
     z1 = tmpZ
     kb = tmpk
@@ -179,6 +199,10 @@ subroutine displacement(xd, yd, zd, ib, jb, kb, dt,imt,jmt,h,ah,av,dxv,dyu,do3d,
     zd = 0.
 endif
 
+! print *,' '
+! print *,'circle'
+! print *,'xd=',xd,' yd=',yd
+
 if(doturb==3) then
     !_______________________________________________________________________________
     ! The diffusion is here set on an ellipse instead of a circle
@@ -259,6 +283,8 @@ subroutine displacement(xd, yd, zd, ib, jb, kb, dt,imt,jmt,h,ah,av,dxv,dyu,do3d,
     xd= xx*dcos(theta)-yy*dsin(theta)
     yd=-xx*dsin(theta)+yy*dcos(theta)
 
+!     print *,'ellipse'
+!     print *,'xd=',xd,' yd=',yd
     ! print *,'theta=',theta,' xx=',xx,' yy=',yy,' grdx=',grdx,' grdy=',grdy
 
     !print *,'theta=',theta*180./pi,elip,grdx/grad

init.py

@@ -65,7 +65,8 @@
 import glob
 from datetime import datetime, timedelta
 from matplotlib.mlab import *
-
+import inout
+import tools
 
 def galveston():
 	'''
@@ -126,37 +127,59 @@ def galveston():
 def test1():
 	'''
 	A drifter test using TXLA model output. 
-	The comparison case for this simulation is 2D with no turbulence/diffusion.
+	The comparison case for this simulation is 2D (do3d=0) 
+	with no turbulence/diffusion (doturb=0).
 	Drifters are started at the surface and run forward
-	for five days from 11/25/09. Compare results with figure in examples/test1.png.
+	for ten days (ndays=10) from 11/25/09 (in date). Compare results with figure in examples/test1.png.
 
 	'''
 
 	# Location of TXLA model output
-	# file and then grid
-	loc = ['http://barataria.tamu.edu:8080/thredds/dodsC/txla_nesting6/ocean_his_0150.nc', \
-			'http://barataria.tamu.edu:8080//thredds/dodsC/txla_nesting6_grid/txla_grd_v4_new.nc']
+	# file and then grid. 
+	# 0150 file goes from (2009, 11, 19, 12, 0) to (2009, 12, 6, 0, 0)
+	# loc = ['http://barataria.tamu.edu:8080/thredds/dodsC/txla_nesting6/ocean_his_0150.nc', \
+	# 		'http://barataria.tamu.edu:8080//thredds/dodsC/txla_nesting6_grid/txla_grd_v4_new.nc']
+	# Location of TXLA model output
+	if 'rainier' in os.uname():
+		loc = '/Users/kthyng/Documents/research/postdoc/' # for model outputs
+	elif 'hafen.tamu.edu' in os.uname():
+		loc = '/home/kthyng/shelf/' # for model outputs
 
 	# Initialize parameters
-	nsteps = 10
-	ndays = 2
+	nsteps = 5
+	ndays = 10 #16
 	ff = 1
 	# Start date
 	date = datetime(2009,11, 25, 0)
+	# date = datetime(2009,11, 20, 0)
+
 	# Time between outputs
 	# Dt = 14400. # in seconds (4 hours), nc.variables['dt'][:] 
 	tseas = 4*3600 # 4 hours between outputs, in seconds, time between model outputs 
-	ah = 100.
+	ah = 5. #100.
 	av = 1.e-5 # m^2/s, or try 5e-6
 
+	# grid = netCDF.Dataset(loc+'grid.nc')
+	# lonr = grid.variables['lon_rho'][:]
+	# latr = grid.variables['lat_rho'][:]
+	grid = inout.readgrid(loc)
+
 	## Input starting locations as real space lon,lat locations
 	# lon0,lat0 = np.meshgrid(-95.498218005315309,23.142258627126882) # [0,0] (SE) corner
 	# lon0,lat0 = np.meshgrid(-97.748582291691989,23.000027311710628) # [-1,0] (SW) corner
 	# lon0,lat0 = np.meshgrid(-87.757124031927574,29.235771320764623) # [0,-1] (NE) corner
 	# lon0,lat0 = np.meshgrid(-88.3634073986196,30.388542615201313) # [-1,-1] (NW) corner
-	lon0,lat0 = np.meshgrid(np.linspace(-94,-93,5),np.linspace(28,29,5)) # grid outside Galveston Bay
-	lon0 = lon0.flatten()
-	lat0 = lat0.flatten()
+	# lon0,lat0 = np.meshgrid(np.linspace(-94,-93,10),np.linspace(28,29,10)) # grid outside Galveston Bay
+	# lon0,lat0 = np.meshgrid(np.linspace(-95,-91,100),np.linspace(28,29,50)) # rectangle outside Galveston
+
+	# lon0,lat0 = np.meshgrid(np.linspace(-98.5,-87.5,1100),np.linspace(22.5,31,980)) # whole domain, 1 km
+	# lon0,lat0 = np.meshgrid(np.linspace(-98.5,-87.5,220),np.linspace(22.5,31,196)) # whole domain, 5 km
+	# FOR TEST1:
+	lon0,lat0 = np.meshgrid(np.linspace(-98.5,-87.5,110),np.linspace(22.5,31,98)) # whole domain, 10 km
+	# lon0,lat0 = np.meshgrid(np.linspace(-98.5,-87.5,21),np.linspace(22.5,31,20)) # whole domain, 50 km
+
+	# Eliminate points that are outside domain or in masked areas
+	lon0,lat0 = tools.check_points(lon0,lat0,grid)
 
 	## Choose method for vertical placement of drifters
 	# Also update makefile accordingly. Choose the twodim flag for isoslice.
@@ -177,15 +200,20 @@ def test1():
 	doturb = 0
 
 	# simulation name, used for saving results into netcdf file
-	name = 'test1'
+	name = 'test1' #'5_5_D5_F'
+
+	# Save these settings to a file
+	np.savez('tracks/' + name + 'header.in',loc=loc,nsteps=nsteps,ndays=ndays, 
+			ff=ff,date=date,tseas=tseas,ah=ah,av=av,lon0=lon0,lat0=lat0,
+			z0=z0,zpar=zpar,do3d=do3d,doturb=doturb,name=name)
 
 	return loc,nsteps,ndays,ff,date,tseas,ah,av,lon0,lat0,z0,zpar,do3d,doturb,name
 
 def test2():
 	'''
 	A drifter test using TXLA model output. 
-	This simulation is 3D with turbulence (doturb=1) added in.
-	Drifters are started at 10 meters below the mean sea level and run backward
+	This simulation is 3D (do3d=1) with turbulence (doturb=1) added in.
+	Drifters are started at 10 meters below the mean sea level and run backward (ff=-1)
 	for five days from 11/25/09. Compare results with figure in examples/test2.png.
 	'''
 
@@ -197,7 +225,7 @@ def test2():
 	# Initialize parameters
 	nsteps = 10
 	ndays = 5
-	ff = 1
+	ff = -1
 	# Start date
 	date = datetime(2009,11, 25, 0)
 	# Time between outputs