Merge pull request #1000 from MJHarrison-GFDL/ALE_sponges_ongrid

Ale sponges ongrid

.testing/_tc4/build_data.py

@@ -0,0 +1,68 @@
+import netCDF4 as nc
+import numpy as np
+
+x=nc.Dataset('ocean_hgrid.nc').variables['x'][1::2,1::2]
+y=nc.Dataset('ocean_hgrid.nc').variables['y'][1::2,1::2]
+zbot=nc.Dataset('topog.nc').variables['depth'][:]
+zbot0=zbot.max()
+
+def t_fc(x,y,z,radius=5.0,tmag=1.0):  # a radially symmetric anomaly in the center of the domain. units are meters and degC
+    ny,nx=x.shape;nz=z.shape[0]
+    x0=x[int(ny/2),int(nx/2)];y0=y[int(ny/2),int(nx/2)]
+    tl=np.zeros((nz,ny,nx))
+    zb=z[-1]
+    if len(z)>1:
+        zd=z/zb
+    else:
+        zd=[0.]
+    for k in np.arange(len(zd)):
+        r=np.sqrt((x-x0)**2.+(y-y0)**2.)
+        tl[k,:]=tl[k,:]+(1.0-np.minimum(r/radius,1.0))*tmag*(1.0-zd[k])
+    return tl
+
+ny,nx = x.shape
+nz=10;z=(np.arange(nz)*zbot0)/nz
+
+temp=t_fc(x,y,z)
+salt=np.zeros(temp.shape)+35.0
+fl=nc.Dataset('temp_salt_ic.nc','w',format='NETCDF3_CLASSIC')
+fl.createDimension('lon',nx)
+fl.createDimension('lat',ny)
+fl.createDimension('depth',nz)
+fl.createDimension('Time',None)
+zv=fl.createVariable('depth','f8',('depth'))
+lonv=fl.createVariable('lon','f8',('lon'))
+latv=fl.createVariable('lat','f8',('lat'))
+timev=fl.createVariable('Time','f8',('Time'))
+timev.calendar='noleap'
+timev.units='days since 0001-01-01 00:00:00.0'
+timev.modulo=' '
+tv=fl.createVariable('ptemp','f8',('Time','depth','lat','lon'),fill_value=-1.e20)
+sv=fl.createVariable('salt','f8',('Time','depth','lat','lon'),fill_value=-1.e20)
+tv[:]=temp[np.newaxis,:]
+sv[:]=salt[np.newaxis,:]
+zv[:]=z
+lonv[:]=x[0,:]
+latv[:]=y[:,0]
+timev[0]=0.
+fl.sync()
+fl.close()
+
+
+# Make Sponge forcing file
+dampTime=20.0 # days
+secDays=8.64e4
+fl=nc.Dataset('sponge.nc','w',format='NETCDF3_CLASSIC')
+fl.createDimension('lon',nx)
+fl.createDimension('lat',ny)
+lonv=fl.createVariable('lon','f8',('lon'))
+latv=fl.createVariable('lat','f8',('lat'))
+spv=fl.createVariable('Idamp','f8',('lat','lon'),fill_value=-1.e20)
+Idamp=np.zeros((ny,nx))
+if dampTime>0.:
+    Idamp=0.0+1.0/(dampTime*secDays)
+spv[:]=Idamp
+lonv[:]=x[0,:]
+latv[:]=y[:,0]
+fl.sync()
+fl.close()

.testing/_tc4/build_grid.py

@@ -0,0 +1,75 @@
+import netCDF4 as nc
+from netCDF4 import stringtochar
+import numpy as np
+
+
+nx=14;ny=10 # grid size
+depth0=100. #uniform depth
+ds=0.01 # grid resolution at the equator in degrees
+Re=6.378e6 # Radius of earth
+
+topo_=np.zeros((ny,nx))+depth0
+f_topo=nc.Dataset('topog.nc','w',format='NETCDF3_CLASSIC')
+ny,nx=topo_.shape
+f_topo.createDimension('ny',ny)
+f_topo.createDimension('nx',nx)
+f_topo.createDimension('ntiles',1)
+f_topo.createVariable('depth','f8',('ny','nx'))
+f_topo.createVariable('h2','f8',('ny','nx'))
+f_topo.variables['depth'][:]=topo_
+f_topo.sync()
+f_topo.close()
+
+x_=np.arange(0,2*nx+1)*ds # units are degrees E
+y_=np.arange(0,2*ny+1)*ds # units are degrees N
+x,y=np.meshgrid(x_,y_)
+
+dx=np.zeros((2*ny+1,2*nx))
+dy=np.zeros((2*ny,2*nx+1))
+rad_deg=np.pi/180.
+dx[:]=rad_deg*Re*(x[:,1:]-x[:,0:-1])*np.cos(0.5*rad_deg*(y[:,0:-1]+y[:,1:]))
+dy[:]=rad_deg*Re*(y[1:,:]-y[0:-1,:])
+
+f_sg=nc.Dataset('ocean_hgrid.nc','w',format='NETCDF3_CLASSIC')
+f_sg.createDimension('ny',ny*2)
+f_sg.createDimension('nx',nx*2)
+f_sg.createDimension('nyp',ny*2+1)
+f_sg.createDimension('nxp',nx*2+1)
+f_sg.createDimension('string',5)
+f_sg.createVariable('y','f8',('nyp','nxp'))
+f_sg.createVariable('x','f8',('nyp','nxp'))
+dyv=f_sg.createVariable('dy','f8',('ny','nxp'))
+dxv=f_sg.createVariable('dx','f8',('nyp','nx'))
+areav=f_sg.createVariable('area','f8',('ny','nx'))
+dxv.units='m'
+dyv.units='m'
+areav.units='m2'
+f_sg.createVariable('angle_dx','f8',('nyp','nxp'))
+f_sg.createVariable('tile','S1',('string'))
+f_sg.variables['y'].units='degrees'
+f_sg.variables['x'].units='degrees'
+f_sg.variables['dy'].units='meters'
+f_sg.variables['dx'].units='meters'
+f_sg.variables['area'].units='m2'
+f_sg.variables['angle_dx'].units='degrees'
+f_sg.variables['y'][:]=y
+f_sg.variables['x'][:]=x
+f_sg.variables['dx'][:]=dx
+f_sg.variables['dy'][:]=dy
+#Compute the area bounded by lines of constant
+#latitude-longitud on a sphere in m2.
+dlon=x_[1:]-x_[:-1]
+dlon=np.tile(dlon[np.newaxis,:],(2*ny,1))
+y1_=y_[:-1]
+y1_=y1_[:,np.newaxis]*rad_deg
+y2_=y_[1:]
+y2_=y2_[:,np.newaxis]*rad_deg
+y1_=np.tile(y1_,(1,2*nx))
+y2_=np.tile(y2_,(1,2*nx))
+area=(rad_deg*Re*Re)*(np.sin(y2_)-np.sin(y1_)) * dlon
+f_sg.variables['area'][:]=area
+f_sg.variables['angle_dx'][:]=0.
+str_=stringtochar(np.array(['tile1'],dtype='S5'))
+f_sg.variables['tile'][:] = str_
+f_sg.sync()
+f_sg.close()

.testing/_tc4/diag_table

@@ -0,0 +1,2 @@
+"MOM test configuration 4"
+1 1 1 0 0 0

.testing/_tc4/input.nml

@@ -0,0 +1,27 @@
+ &MOM_input_nml
+         output_directory = './',
+         input_filename = 'n'
+         restart_input_dir = 'INPUT/',
+         restart_output_dir = 'RESTART/',
+         parameter_filename = 'MOM_input',
+                              'MOM_override' /
+
+ &diag_manager_nml
+         flush_nc_files = .true.
+ /
+
+ &fms_nml
+         domains_stack_size = 710000,
+         stack_size = 0 /
+
+ &ocean_domains_nml
+ /
+
+ &ocean_solo_nml
+            months = 0
+            date_init = 1,1,1,0,0,0
+            hours = 0
+            minutes = 0
+            seconds = 0
+            calendar = 'julian' /
+

src/framework/MOM_horizontal_regridding.F90

@@ -734,16 +734,15 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
 
   if (z_edges_in(kd+1)<max_depth) z_edges_in(kd+1)=max_depth
 
-  if (is_root_pe()) &
-    call time_interp_external(fms_id, Time, data_in, verbose=.true.)
-
   !  roundoff = 3.0*EPSILON(missing_value)
   roundoff = 1.e-4
 
   if (.not.spongeDataOngrid) then
-  ! loop through each data level and interpolate to model grid.
-  ! after interpolating, fill in points which will be needed
-  ! to define the layers
+    if (is_root_pe()) &
+      call time_interp_external(fms_id, Time, data_in, verbose=.true.)
+    ! loop through each data level and interpolate to model grid.
+    ! after interpolating, fill in points which will be needed
+    ! to define the layers
     do k=1,kd
       write(laynum,'(I8)') k ; laynum = adjustl(laynum)
       if (is_root_pe()) then
@@ -860,6 +859,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
 
     enddo ! kd
   else
+      call time_interp_external(fms_id, Time, data_in, verbose=.true.)
       do k=1,kd
         do j=js,je
           do i=is,ie

src/parameterizations/vertical/MOM_ALE_sponge.F90

@@ -902,7 +902,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
       sp_val(:,:,:)=0.0
       mask_z(:,:,:)=0.0
 
-      call horiz_interp_and_extrap_tracer(CS%Ref_val(CS%fldno)%id,Time, 1.0,G,sp_val,mask_z,z_in,z_edges_in, &
+      call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id,Time, 1.0,G,sp_val,mask_z,z_in,z_edges_in, &
                       missing_value,.true., .false.,.false., m_to_Z=US%m_to_Z,spongeOnGrid=CS%SpongeDataOngrid)
 
 !     call pass_var(sp_val,G%Domain)