modifying the benchmark accordingly

benchmarks/titan/ocean/functions.py

@@ -35,4 +35,4 @@ def tideDynColat(lonlat, t):
 
   dObldColat = obl*p_const/R*np.cos(2.*colat)*np.cos(lon)*np.cos(omega*t) # 1/R *dphi_ob/dcolat
 
-  return -(dObldColat + dEccdColat)*-1# *-1 because lat and colat are anti parallel
+  return -(dObldColat + dEccdColat)

benchmarks/titan/ocean/prepro.py

@@ -46,22 +46,28 @@ if "coriolis" in do:
     slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'coriolis.nc', 'coriolis', param.data_dir+'coriolis')  
     
 if "forcing" in do:
-  print('Preprocessing Forcing')
+ print('Preprocessing Forcing')
   lonlat = lonlat_global.coordinates
   nPt = len(lonlat)
-  
+  Flon = np.empty(nPt)
+  Fcolat = np.empty(nPt)
+  #Fr = np.zeros(nPt)
+         
   Fx = np.empty((param.nt_prepro,nPt))
   Fy = np.empty((param.nt_prepro,nPt))
   Fz = np.empty((param.nt_prepro,nPt))
+
   time_vec = np.empty((param.nt_prepro))
-    
+
   t = param.initial_time
   for i in range (param.nt_prepro):
     for j in range(nPt):
-      Flon= functions.tideDynLon(lonlat[j,:], t)
-      Fcolat = functions.tideDynColat(lonlat[j,:], t)
+      Flon[j] = functions.tideDynLon(lonlat[j,:], t)
+      Fcolat[j] = functions.tideDynColat(lonlat[j,:], t)     
       Fr = 0
-      Fx[i,j], Fy[i,j], Fz[i,j] = lonlat_global.rotate3D(Flon, Fcolat, Fr,lonlat[j,:])
+      #Fx[i,j], Fy[i,j], Fz[i,j] = lonlat_global.rotate3D(Flon[j], Fcolat[j], Fr,lonlat[j,:])      
+    (Fx[i,:], Fy[i,:], Fz[i,:]) = lonlat_global.rotateSphere(Flon, Fcolat, 0)
+                                                
     time_vec[i] = t
     turn = t // param.T
     turnRest = float(t % param.T)
@@ -69,12 +75,13 @@ if "forcing" in do:
     print('|ITER|',i,'|TURN|',turn,'|DAYS|',days,'|DT|',param.dt, '|t|', t)
     t = t + param.dt
     i=i+1
-
   time = slimPre.Time(time_vector = time_vec, periodic = True)
-  slimPre.write_file(param.data_dir+'forcing.nc', region=region_global, time=time, data=[('Fx', Fx), ('Fy', Fy), ('Fz', Fz)])
+  slimPre.write_file(param.data_dir+'forcing_'+param.bathy_acc+'.nc', region=region_global, time=time, data=[('Fx', Fx), ('Fy', Fy), ('Fz', Fz)])
   if check_result:
-    slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing.nc', 'Fx', param.data_dir+'Fx',50*param.dt) 
-    slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing.nc', 'Fy', param.data_dir+'Fy',50*param.dt) 
+    for i in range (param.nt_prepro):
+      slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing_'+param.bathy_acc+'.nc', 'Fx', param.data_dir+'Fx-%03d' %i,i*param.dt) 
+      slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing_'+param.bathy_acc+'.nc', 'Fy', param.data_dir+'Fy-%03d' %i,i*param.dt) 
+      slimPre.netcdf_to_msh(mesh_file_name, param.data_dir+'forcing_'+param.bathy_acc+'.nc', 'Fz', param.data_dir+'Fz-%03d' %i,i*param.dt)
     
 if "gravity" in do:
   lonlat = lonlat_global.coordinates