Skip to content

GitLab

Commit 868014cc authored by ashbre's avatar ashbre
Browse files

Adding structure to namelists and fortran files and other small files

parent ed6b5111
Hide whitespace changes
Inline Side-by-side
Showing with 5675 additions and 0 deletions
NAMELISTS_AND_FORTRAN_FILES/INITIAL_CONDITION/initcd_vosaline.namelist 0 → 100755
View file @ 868014cc
!! -------------------
!! Namelist for SOSIE
!! -------------------
!!
!!
!! ************************************************************
!! &ninput => info about field to interpolate
!! and source grid to interpolate from
!! ************************************************************
!!
!! ivect : vector correction control on treated field [integer]
!! ivect = 0 : input field is not a component of a vector
!! or the target grid is regular (lregout = T)
!! * if non-regular distorted target grids (like ORCAX):
!! ivect = 1 : input field is a zonal (X) component of a vector
!! ivect = 2 : input field is a meridional (Y) component of a vector
!!
!! lregin : is the source grid regular? [logical]
!! (ie : are input longitude and latitude 1D?)
!!
!! cf_in : file containing the input field to be interpolated [char]
!! cv_in : name of treated variable (in input field file) [char]
!!
!! cv_t_in : name of time record variable in the input file [char]
!! or 'missing_rec' if no time record is present in the input file
!!
!! jt1 : first time record to be interpolated
!! jt2 : last time record to be interpolated
!! => set jt1 and jt2 to 0 if you want to skip this option
!! and interpolate the nt time steps of the current field
!!
!! jplev : level to treat if your file is 3D (spatial), has no influence if
!! your file is 2D in space !
!! ------------------------------------------------------------------
!! jplev > 0 = level to treat (ex : jplev = 1 will interpolate only
!! surface field corresponding to the 1st level )
!! ------------------------------------------------------------------
!! jplev = 0 : perform 3D interpolation (if input file is 3D) !!! |
!! ------------------------------------------------------------------
!! jplev = -1 : overrides good sense and forces sosie to understand that
!! your field to interpolate is 2D with a time record
!! (usually the case if the time record dimension in your
!! input file is not declared as UNLIMITED => bad! )
!! => so SOSIE doesn't mistake this time record with a depth!
!! -------------------------------------------------------------------
!!
!! cf_x_in : file containing the input grid (usually = cf_in) [char]
!! cv_lon_in : name of longitude in the input grid file [char]
!! cv_lat_in : name of latitude in the input grid file [char]
!!
!! cf_lsm_in : (only relevant if ldrown==.true.)
!! file containing the input land-sea mask [char]
!! Alternatively:
!! * specify " cf_lsm_in = 'missing_value' " if a 'missing_value' netcdf
!! attribute defines the mask on the input data field
!! * specify " cf_lsm_in = 'nan' " if mask is defined with NaN values
!! * specify " cf_lsm_in = 'value' if you want land regions to be defined
!! where field 'cv_in' is strictly equal to the numeric value read into 'cv_lsm_in'
!! * specify " cf_lsm_in = 'val+' if you want land regions to be defined
!! where field 'cv_in' is larger or equal to the numeric value read into 'cv_lsm_in'
!! * specify " cf_lsm_in = 'val-' if you want land regions to be defined
!! where field 'cv_in' is smaller or equal to the numeric value read into 'cv_lsm_in'
!! Ex: you want all points where your field is <= 0 to become land mask,
!! then specify: cf_lsm_in = 'val-' and cv_lsm_in = '0.00001'
!!
!! cv_lsm_in : (only relevant if ldrown==.true.)
!! name of land-sea mask variable [char]
!! or if cf_lsm_in = 'missing_value'--> '')
!! by default ocean is flagged with value 1
!! and continents are flagged with value 0
!! Alternatively:
!! a string of numeric value when cf_lsm_in is 'value', 'val-', or 'val+'
!!
!! ldrown : whether we call DROWN land filling procedure [logical]
!! => will propagate/extrapolate sea values (defined where lsm==1)
!! of field cv_in ONTO continents (defined WHERE lsm==0) to avoid
!! interpolation problems, such as continental values that contaminate
!! sea values during interpolation
!!
!! ewper : east-west periodicity on the input file/grid [integer]
!! = -1 --> no periodicity
!! >= 0 --> periodicity with overlap of ewper points
!!
!! vmax : upper bound not to exceed for treated variable [real]
!! vmin : lower bound not to exceed for treated variable [real]
!!
!! ismooth : if ismooth > 0 the field to be interpolated will be smoothed
!! prior to interpolation. By applying ismooth times a type of
!! closest neighboors boxcar smoothing algorithm
!! (check "SMOOTH" of mod_drown.f90)
!! => this is usefull to avoid sub-sampling when your target
!! grid is much coarser than your source grid
!! (i.e. when interpolating from high-res to low-res)
!! => start with a multiple of 10, typically 20, and adjust depending
!! on the result
!!
!!--------------------------------------------------------------------------
!!
&ninput
ivect = 0
lregin = T
cf_in = './initial_condition.nc'
cv_in = 'so'
cv_t_in = 'time'
jt1 = 0
jt2 = 0
jplev = 0
cf_x_in = './initial_condition.nc'
cv_lon_in = 'longitude'
cv_lat_in = 'latitude'
cf_lsm_in = 'missing_value'
cv_lsm_in = 'mask'
ldrown = T
ewper = -1
vmax = 1.E6
vmin = -1.E6
ismooth = 0
/
!!
!!
!!
!!
!! ***********************************************************
!! &n3d => info about source and target vertical levels/grids
!! ONLY IF 3D INTERPOLATION ( jplev = 0 in &ninput)
!! ***********************************************************
!!
!! Only mind if you do want to perform a 3D (spatial) interpolation
!!
!! Mind only if you do want to perform a 3D interpolation !
!! First, make sure that jplev is set to 0 !
!!
!! cf_z_in : file containing the input depth vector (associates a depth to a
!! given level). In most cases should be the same file than cf_x_in.
!! cv_z_in : name of the variable for the input depth vector
!!
!! cf_z_out : file containing the output depth vector (associates a depth to a
!! given level). In most cases should be the same file than cf_x_in.
!! cv_z_out : name of the variable for the output depth vector in file 'cf_z_out'
!! cv_z_out_name: name you wish to give to depth variable in file to be created...
!!
!! ctype_z_in : type of coordinates in input file (currently available z/sigma)
!! ctype_z_out : type of coordinates in output file (currently available z/sigma)
!!
!! These are to be set ONLY if ctype_z_in = 'sigma'
!! cf_bathy_in : file containing the bathymetry on input grid (usually ROMS grid file)
!! cv_bathy_in : name of bathymetry variable (usually h)
!! ssig_in : structure with ROMS s-coordinates parameters on input grid
!! Vtransform | Vstretching | Nlevels | theta_s | theta_b | Tcline | hmin
!!
!! These are to be set ONLY if ctype_z_out = 'sigma'
!! cf_bathy_out : file containing the bathymetry on output grid (usually ROMS grid file)
!! cv_bathy_out : name of bathymetry variable (usually h)
!! ssig_out : structure with ROMS s-coordinates parameters on output grid (see above)
!!
&n3d
cf_z_in = 'initial_condition.nc'
cv_z_in = 'depth'
cf_z_out = 'domain_cfg.nc'
cv_z_out = 'nav_lev'
cv_z_out_name = 'gdept'
ctype_z_in = 'z'
ctype_z_out = 'z'
/
!!
!!
!!
!!
!!
!! *****************************************************************
!! &nhtarget => info about horizontal target grid to interpolate to
!! *****************************************************************
!!
!! lregout : is the target grid regular ? [logical]
!! (ie : are output longitude and latitude 1D?)
!!
!! cf_x_out : file containing the target grid [char]
!! cv_lon_out : name of longitude variable [char]
!! cv_lat_out : name of latitude variable [char]
!!
!! TRICK: for interpolating onto a global regular spherical grid
!! ------ with a resolution of dx deg. of longitude and dy deg. of latitude
!! * cf_x_out = 'spheric' ! tells SOSIE to build a spherical output grid
!! * cv_lon_out = '1.0' ! your dx, here 1.0 deg.
!! * cv_lat_out = '1.0' ! your dy, here 1.0 deg.
!!
!!
!! cf_lsm_out : file containing output land-sea mask [char]
!! MUST BE 3D for 3D interpolation!
!! or specify 'missing_value' if a 'missing_value' netcdf
!! attribute defines the mask on a field 'X' in file 'cf_x_out'
!! (not needed if "lmout = .FALSE." --> '')
!!
!! cv_lsm_out : name of land-sea mask variable in 'cf_lsm_out' [char]
!! or name of field 'X' in 'cf_x_out' if you specified
!! cf_lsm_out = 'missing_value'
!! (not needed if "lmout = .FALSE." --> '')
!!
!! lmout : whether to mask the interpolated field on the output file [logical]
!! if lmout is set to .FALSE. and cf_lsm_out is different than '' the output
!! field will be drowned using the mask defined by cf_lsm_out (and cv_lsm_out)
!!
!! rmaskvalue : missing value given to output field (for continents) [logical]
!!
!! lct : whether to control or not time variable [logical]
!! TRUE -> specify time array with starting time 't0' and step 't_stp'
!! usefull if you do not have a "time" variable in your input netcdf file !
!! FALSE -> same time array as in input file is used
!! t0 : time to start (if lct is set to .TRUE.) [real]
!! t_stp : time step (if lct is set to .TRUE.) [real]
!!
!! ewper_out : east-west periodicity on the output file/grid [integer]
!! = -1 --> no periodicity
!! >= 0 --> periodicity with overlap of ewper points
!!
!!
&nhtarget
lregout = T
cf_x_out = 'initial_condition.nc'
cv_lon_out = 'longitude'
cv_lat_out = 'latitude'
cf_lsm_out = ''
cv_lsm_out = ''
lmout = F
!rmaskvalue = -9999
lct = F
t0 = 0.
t_stp = 0.
ewper_out = -1
/
!!
!!
!!
!!
!! *****************************************************************
!! &noutput => info on the (horizontal) interpolation method to use
!! and the netcdf file to generate
!! *****************************************************************
!!
!! This mostly deals with how the output file to be created is going to look like!
!!
!! cmethod : the 2D interpolation method to be used
!!
!! * use 'akima' if your input domain is regular (non-distorted grid)
!!
!! * use 'bilin' otherwise (bilinear 2D interpolation)
!!
!! * use 'no_xy' to only perform vertical interpolation, i.e. interpolate a
!! a 3D field given on ni*nj and nk_in levels to the same ni*nj 2D domain
!! but on nk_out levels!
!! => for example interpolates a 3D field from grid ORCAX.L46 to ORCAX.L75
!!
!! *** Into the netcdf file to be created : ***
!! cv_t_out : name of time record vector in the output file [char]
!! => set to cv_t_out='' if no time dimension
!! cv_out : name for treated variable in the output file [char]
!! cu_out : if not = '': then change the unit of treated variable units [char]
!! cln_out : if not = '': then change the long name treated variable [char]
!! cd_out : directory to create output file to [char]
!!
!! *** Naming of the output file : ***
!! csource : short string to describe the origin grid [char]
!! ctarget : short string to describe the target grid [char]
!! cextra : short extra indication about the file [char]
!!
&noutput
cmethod = 'bilin'
cv_t_out = 'time_counter'
cv_out = 'vosaline'
cu_out = 'PSU'
cln_out = 'Salinity'
cd_out = '.'
csource = 'COPERNICUS'
ctarget = 'INDIAN'
cextra = '2016'
/
!!
NAMELISTS_AND_FORTRAN_FILES/INITIAL_CONDITION/initcd_votemper.namelist 0 → 100755
View file @ 868014cc
!! -------------------
!! Namelist for SOSIE
!! -------------------
!!
!!
!! ************************************************************
!! &ninput => info about field to interpolate
!! and source grid to interpolate from
!! ************************************************************
!!
!! ivect : vector correction control on treated field [integer]
!! ivect = 0 : input field is not a component of a vector
!! or the target grid is regular (lregout = T)
!! * if non-regular distorted target grids (like ORCAX):
!! ivect = 1 : input field is a zonal (X) component of a vector
!! ivect = 2 : input field is a meridional (Y) component of a vector
!!
!! lregin : is the source grid regular? [logical]
!! (ie : are input longitude and latitude 1D?)
!!
!! cf_in : file containing the input field to be interpolated [char]
!! cv_in : name of treated variable (in input field file) [char]
!!
!! cv_t_in : name of time record variable in the input file [char]
!! or 'missing_rec' if no time record is present in the input file
!!
!! jt1 : first time record to be interpolated
!! jt2 : last time record to be interpolated
!! => set jt1 and jt2 to 0 if you want to skip this option
!! and interpolate the nt time steps of the current field
!!
!! jplev : level to treat if your file is 3D (spatial), has no influence if
!! your file is 2D in space !
!! ------------------------------------------------------------------
!! jplev > 0 = level to treat (ex : jplev = 1 will interpolate only
!! surface field corresponding to the 1st level )
!! ------------------------------------------------------------------
!! jplev = 0 : perform 3D interpolation (if input file is 3D) !!! |
!! ------------------------------------------------------------------
!! jplev = -1 : overrides good sense and forces sosie to understand that
!! your field to interpolate is 2D with a time record
!! (usually the case if the time record dimension in your
!! input file is not declared as UNLIMITED => bad! )
!! => so SOSIE doesn't mistake this time record with a depth!
!! -------------------------------------------------------------------
!!
!! cf_x_in : file containing the input grid (usually = cf_in) [char]
!! cv_lon_in : name of longitude in the input grid file [char]
!! cv_lat_in : name of latitude in the input grid file [char]
!!
!! cf_lsm_in : (only relevant if ldrown==.true.)
!! file containing the input land-sea mask [char]
!! Alternatively:
!! * specify " cf_lsm_in = 'missing_value' " if a 'missing_value' netcdf
!! attribute defines the mask on the input data field
!! * specify " cf_lsm_in = 'nan' " if mask is defined with NaN values
!! * specify " cf_lsm_in = 'value' if you want land regions to be defined
!! where field 'cv_in' is strictly equal to the numeric value read into 'cv_lsm_in'
!! * specify " cf_lsm_in = 'val+' if you want land regions to be defined
!! where field 'cv_in' is larger or equal to the numeric value read into 'cv_lsm_in'
!! * specify " cf_lsm_in = 'val-' if you want land regions to be defined
!! where field 'cv_in' is smaller or equal to the numeric value read into 'cv_lsm_in'
!! Ex: you want all points where your field is <= 0 to become land mask,
!! then specify: cf_lsm_in = 'val-' and cv_lsm_in = '0.00001'
!!
!! cv_lsm_in : (only relevant if ldrown==.true.)
!! name of land-sea mask variable [char]
!! or if cf_lsm_in = 'missing_value'--> '')
!! by default ocean is flagged with value 1
!! and continents are flagged with value 0
!! Alternatively:
!! a string of numeric value when cf_lsm_in is 'value', 'val-', or 'val+'
!!
!! ldrown : whether we call DROWN land filling procedure [logical]
!! => will propagate/extrapolate sea values (defined where lsm==1)
!! of field cv_in ONTO continents (defined WHERE lsm==0) to avoid
!! interpolation problems, such as continental values that contaminate
!! sea values during interpolation
!!
!! ewper : east-west periodicity on the input file/grid [integer]
!! = -1 --> no periodicity
!! >= 0 --> periodicity with overlap of ewper points
!!
!! vmax : upper bound not to exceed for treated variable [real]
!! vmin : lower bound not to exceed for treated variable [real]
!!
!! ismooth : if ismooth > 0 the field to be interpolated will be smoothed
!! prior to interpolation. By applying ismooth times a type of
!! closest neighboors boxcar smoothing algorithm
!! (check "SMOOTH" of mod_drown.f90)
!! => this is usefull to avoid sub-sampling when your target
!! grid is much coarser than your source grid
!! (i.e. when interpolating from high-res to low-res)
!! => start with a multiple of 10, typically 20, and adjust depending
!! on the result
!!
!!--------------------------------------------------------------------------
!!
&ninput
ivect = 0
lregin = T
cf_in = './initial_condition.nc'
cv_in = 'thetao'
cv_t_in = 'time'
jt1 = 0
jt2 = 0
jplev = 0
cf_x_in = './initial_condition.nc'
cv_lon_in = 'longitude'
cv_lat_in = 'latitude'
cf_lsm_in = 'missing_value'
cv_lsm_in = 'mask'
ldrown = T
ewper = -1
vmax = 1.E6
vmin = -1.E6
ismooth = 0
/
!!
!!
!!
!!
!! ***********************************************************
!! &n3d => info about source and target vertical levels/grids
!! ONLY IF 3D INTERPOLATION ( jplev = 0 in &ninput)
!! ***********************************************************
!!
!! Only mind if you do want to perform a 3D (spatial) interpolation
!!
!! Mind only if you do want to perform a 3D interpolation !
!! First, make sure that jplev is set to 0 !
!!
!! cf_z_in : file containing the input depth vector (associates a depth to a
!! given level). In most cases should be the same file than cf_x_in.
!! cv_z_in : name of the variable for the input depth vector
!!
!! cf_z_out : file containing the output depth vector (associates a depth to a
!! given level). In most cases should be the same file than cf_x_in.
!! cv_z_out : name of the variable for the output depth vector in file 'cf_z_out'
!! cv_z_out_name: name you wish to give to depth variable in file to be created...
!!
!! ctype_z_in : type of coordinates in input file (currently available z/sigma)
!! ctype_z_out : type of coordinates in output file (currently available z/sigma)
!!
!! These are to be set ONLY if ctype_z_in = 'sigma'
!! cf_bathy_in : file containing the bathymetry on input grid (usually ROMS grid file)
!! cv_bathy_in : name of bathymetry variable (usually h)
!! ssig_in : structure with ROMS s-coordinates parameters on input grid
!! Vtransform | Vstretching | Nlevels | theta_s | theta_b | Tcline | hmin
!!
!! These are to be set ONLY if ctype_z_out = 'sigma'
!! cf_bathy_out : file containing the bathymetry on output grid (usually ROMS grid file)
!! cv_bathy_out : name of bathymetry variable (usually h)
!! ssig_out : structure with ROMS s-coordinates parameters on output grid (see above)
!!
&n3d
cf_z_in = 'initial_condition.nc'
cv_z_in = 'depth'
cf_z_out = 'domain_cfg.nc'
cv_z_out = 'nav_lev'
cv_z_out_name = 'gdept'
ctype_z_in = 'z'
ctype_z_out = 'z'
/
!!
!!
!!
!!
!!
!! *****************************************************************
!! &nhtarget => info about horizontal target grid to interpolate to
!! *****************************************************************
!!
!! lregout : is the target grid regular ? [logical]
!! (ie : are output longitude and latitude 1D?)
!!
!! cf_x_out : file containing the target grid [char]
!! cv_lon_out : name of longitude variable [char]
!! cv_lat_out : name of latitude variable [char]
!!
!! TRICK: for interpolating onto a global regular spherical grid
!! ------ with a resolution of dx deg. of longitude and dy deg. of latitude
!! * cf_x_out = 'spheric' ! tells SOSIE to build a spherical output grid
!! * cv_lon_out = '1.0' ! your dx, here 1.0 deg.
!! * cv_lat_out = '1.0' ! your dy, here 1.0 deg.
!!
!!
!! cf_lsm_out : file containing output land-sea mask [char]
!! MUST BE 3D for 3D interpolation!
!! or specify 'missing_value' if a 'missing_value' netcdf
!! attribute defines the mask on a field 'X' in file 'cf_x_out'
!! (not needed if "lmout = .FALSE." --> '')
!!
!! cv_lsm_out : name of land-sea mask variable in 'cf_lsm_out' [char]
!! or name of field 'X' in 'cf_x_out' if you specified
!! cf_lsm_out = 'missing_value'
!! (not needed if "lmout = .FALSE." --> '')
!!
!! lmout : whether to mask the interpolated field on the output file [logical]
!! if lmout is set to .FALSE. and cf_lsm_out is different than '' the output
!! field will be drowned using the mask defined by cf_lsm_out (and cv_lsm_out)
!!
!! rmaskvalue : missing value given to output field (for continents) [logical]
!!
!! lct : whether to control or not time variable [logical]
!! TRUE -> specify time array with starting time 't0' and step 't_stp'
!! usefull if you do not have a "time" variable in your input netcdf file !
!! FALSE -> same time array as in input file is used
!! t0 : time to start (if lct is set to .TRUE.) [real]
!! t_stp : time step (if lct is set to .TRUE.) [real]
!!
!! ewper_out : east-west periodicity on the output file/grid [integer]
!! = -1 --> no periodicity
!! >= 0 --> periodicity with overlap of ewper points
!!
!!
&nhtarget
lregout = T
cf_x_out = 'initial_condition.nc'
cv_lon_out = 'longitude'
cv_lat_out = 'latitude'
cf_lsm_out = ''
cv_lsm_out = ''
lmout = F
!rmaskvalue = -9999
lct = F
t0 = 0.
t_stp = 0.
ewper_out = -1
/
!!
!!
!!
!!
!! *****************************************************************
!! &noutput => info on the (horizontal) interpolation method to use
!! and the netcdf file to generate
!! *****************************************************************
!!
!! This mostly deals with how the output file to be created is going to look like!
!!
!! cmethod : the 2D interpolation method to be used
!!
!! * use 'akima' if your input domain is regular (non-distorted grid)
!!
!! * use 'bilin' otherwise (bilinear 2D interpolation)
!!
!! * use 'no_xy' to only perform vertical interpolation, i.e. interpolate a
!! a 3D field given on ni*nj and nk_in levels to the same ni*nj 2D domain
!! but on nk_out levels!
!! => for example interpolates a 3D field from grid ORCAX.L46 to ORCAX.L75
!!
!! *** Into the netcdf file to be created : ***
!! cv_t_out : name of time record vector in the output file [char]
!! => set to cv_t_out='' if no time dimension
!! cv_out : name for treated variable in the output file [char]
!! cu_out : if not = '': then change the unit of treated variable units [char]
!! cln_out : if not = '': then change the long name treated variable [char]
!! cd_out : directory to create output file to [char]
!!
!! *** Naming of the output file : ***
!! csource : short string to describe the origin grid [char]
!! ctarget : short string to describe the target grid [char]
!! cextra : short extra indication about the file [char]
!!
&noutput
cmethod = 'bilin'
cv_t_out = 'time_counter'
cv_out = 'votemper'
cu_out = 'C'
cln_out = 'Temperature'
cd_out = '.'
csource = 'COPERNICUS'
ctarget = 'INDIAN'
cextra = '2016'
/
!!
NAMELISTS_AND_FORTRAN_FILES/INITIAL_CONDITION/namelist_reshape_bilin_initcd_vosaline 0 → 100755
View file @ 868014cc
&comments
-----------------------------------------------------------------------------------
- grid_inputs holds parameters for the scripgrid routine which reformats information
about the input grids
- scripgrid always needs a coordinates.nc file in the
current directory and creates the remapped grid file correspondingly
- it uses the following namelist block to determine its actions
method: only 'regular' is yet implemented, this assumes a cartesian grid
input_lon: name of longitude variable in the input_file
input_lat: name of latitude variable in the input_file
nemo_lon: name of longitude variable in the coordinates.nc
nemo_lat: name of latitude variable in the coordinates.nc
/
&grid_inputs
input_file = 'vosaline_COPERNICUS-INDIAN_2016.nc'
nemo_file = 'coordinates.nc'
datagrid_file = 'remap_data_grid_R12.nc'
nemogrid_file = 'remap_nemo_grid_R12.nc'
method = 'regular'
input_lon = 'longitude'
input_lat = 'latitude'
nemo_lon = 'glamt'
nemo_lat = 'gphit'
nemo_mask = 'none'
nemo_mask_value = 0
input_mask = 'none'
input_mask_value = 0
/
&comments
-----------------------------------------------------------------------------------
- remap_inputs holds parameters for the scrip routine which calculates the weights
needed to convert between two grids
- two remap grid files are required as output by scripgrid
- num_maps is either 1 or 2 depending on whether the reverse transformation is required
- one or two interp_file names are then supplied; these hold the weights to convert
one grid to another
- the map_name variable is just descriptive
- map_method can be 'bilinear' 'conservative' or 'bicubic' (the latter untested)
- normalize_opt should usually be 'frac' or else the user needs to do this scaling
manually (this seems to the case for fractional ice cover)
- restrict_type should be 'latitude' or 'latlon' in which case num_srch_bins only are
used in one or two directions
- use_grid_area fields override the scrip calculation of area in case the model gets
slightly different answers, but the area needs to be supplied in the input files
- output_opt may be supplied and set to either 'scrip' or 'ncar-csm'
/
&remap_inputs
num_maps = 1
grid1_file = 'remap_data_grid_R12.nc'
grid2_file = 'remap_nemo_grid_R12.nc'
interp_file1 = 'data_nemo_bilin_R12.nc'
interp_file2 = 'nemo_data_bilin_R12.nc'
map1_name = 'R12 to nemo bilin Mapping'
map2_name = 'nemo to R12 bilin Mapping'
map_method = 'bilinear'
normalize_opt = 'frac'
output_opt = 'scrip'
restrict_type = 'latitude'
num_srch_bins = 90
luse_grid1_area = .false.
luse_grid2_area = .false.
/
&interp_inputs
input_file = "vosaline_COPERNICUS-INDIAN_2016.nc"
interp_file = "data_nemo_bilin_R12.nc"
input_name = "vosaline"
input_start = 1,1,1,1
input_stride = 1,1,1,1
input_stop = 0,0,0,0
input_vars = "gdept","time_counter"
/
&interp_outputs
output_file = "initcd_vosaline.nc"
output_mode = "create"
output_dims = 'x', 'y', 'z', 'time_counter'
output_scaling = "vosaline|1.0"
output_name = 'vosaline'
output_lon = 'x'
output_lat = 'y'
output_vars = "gdept","time_counter"
&comments
-----------------------------------------------------------------------------------
- shape_inputs holds parameters for the scripshape routine which rearranges the weights
into the form needed by the nemo on the fly interpolation code.
/
&shape_inputs
interp_file = 'data_nemo_bilin_R12.nc'
output_file = 'weights_bilinear_R12.nc'
ew_wrap = -1
/
NAMELISTS_AND_FORTRAN_FILES/INITIAL_CONDITION/namelist_reshape_bilin_initcd_votemper 0 → 100755
View file @ 868014cc
&comments
-----------------------------------------------------------------------------------
- grid_inputs holds parameters for the scripgrid routine which reformats information
about the input grids
- scripgrid always needs a coordinates.nc file in the
current directory and creates the remapped grid file correspondingly
- it uses the following namelist block to determine its actions
method: only 'regular' is yet implemented, this assumes a cartesian grid
input_lon: name of longitude variable in the input_file
input_lat: name of latitude variable in the input_file
nemo_lon: name of longitude variable in the coordinates.nc
nemo_lat: name of latitude variable in the coordinates.nc
/
&grid_inputs
input_file = 'votemper_COPERNICUS-INDIAN_2016.nc'
nemo_file = 'coordinates.nc'
datagrid_file = 'remap_data_grid_R12.nc'
nemogrid_file = 'remap_nemo_grid_R12.nc'
method = 'regular'
input_lon = 'longitude'
input_lat = 'latitude'
nemo_lon = 'glamt'
nemo_lat = 'gphit'
nemo_mask = 'none'
nemo_mask_value = 0
input_mask = 'none'
input_mask_value = 0
/
&comments
-----------------------------------------------------------------------------------
- remap_inputs holds parameters for the scrip routine which calculates the weights
needed to convert between two grids
- two remap grid files are required as output by scripgrid
- num_maps is either 1 or 2 depending on whether the reverse transformation is required
- one or two interp_file names are then supplied; these hold the weights to convert
one grid to another
- the map_name variable is just descriptive
- map_method can be 'bilinear' 'conservative' or 'bicubic' (the latter untested)
- normalize_opt should usually be 'frac' or else the user needs to do this scaling
manually (this seems to the case for fractional ice cover)
- restrict_type should be 'latitude' or 'latlon' in which case num_srch_bins only are
used in one or two directions
- use_grid_area fields override the scrip calculation of area in case the model gets
slightly different answers, but the area needs to be supplied in the input files
- output_opt may be supplied and set to either 'scrip' or 'ncar-csm'
/
&remap_inputs
num_maps = 1
grid1_file = 'remap_data_grid_R12.nc'
grid2_file = 'remap_nemo_grid_R12.nc'
interp_file1 = 'data_nemo_bilin_R12.nc'
interp_file2 = 'nemo_data_bilin_R12.nc'
map1_name = 'R12 to nemo bilin Mapping'
map2_name = 'nemo to R12 bilin Mapping'
map_method = 'bilinear'
normalize_opt = 'frac'
output_opt = 'scrip'
restrict_type = 'latitude'
num_srch_bins = 90
luse_grid1_area = .false.
luse_grid2_area = .false.
/
&interp_inputs
input_file = "votemper_COPERNICUS-INDIAN_2016.nc"
interp_file = "data_nemo_bilin_R12.nc"
input_name = "votemper"
input_start = 1,1,1,1
input_stride = 1,1,1,1
input_stop = 0,0,0,0
input_vars = "gdept","time_counter"
/
&interp_outputs
output_file = "initcd_votemper.nc"
output_mode = "create"
output_dims = 'x', 'y', 'z', 'time_counter'
output_scaling = "votemper|1.0"
output_name = 'votemper'
output_lon = 'x'
output_lat = 'y'
output_vars = "gdept","time_counter"
&comments
-----------------------------------------------------------------------------------
- shape_inputs holds parameters for the scripshape routine which rearranges the weights
into the form needed by the nemo on the fly interpolation code.
/
&shape_inputs
interp_file = 'data_nemo_bilin_R12.nc'
output_file = 'weights_bilinear_R12.nc'
ew_wrap = -1
/
NAMELISTS_AND_FORTRAN_FILES/INITIAL_CONDITION/sosie.x 0 → 100755
View file @ 868014cc
File added
NAMELISTS_AND_FORTRAN_FILES/TIDES/inputs_dst.ncml 0 → 100755
View file @ 868014cc
<ns0:netcdf xmlns:ns0="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2" title="NEMO aggregation">
<ns0:aggregation type="union">
<ns0:netcdf location="file:domain_cfg.nc">
<ns0:variable name="mbathy" orgName="top_level" />
<ns0:variable name="e3u" orgName="e3u_0" />
<ns0:variable name="e3v" orgName="e3v_0" />
</ns0:netcdf>
</ns0:aggregation>
</ns0:netcdf>
NAMELISTS_AND_FORTRAN_FILES/TIDES/inputs_src.ncml 0 → 100755
View file @ 868014cc
<ns0:netcdf xmlns:ns0="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2" title="NEMO aggregation">
<ns0:aggregation type="union">
<ns0:netcdf>
<ns0:aggregation dimName="time_counter" name="temperature" type="joinExisting">
<ns0:netcdf location="http://gws-access.ceda.ac.uk/public/nemo/runs/ORCA0083-N01/means/1979/ORCA0083-N01_19791206d05T.nc" />
</ns0:aggregation>
</ns0:netcdf>
<ns0:netcdf>
<ns0:aggregation dimName="time_counter" name="salinity" type="joinExisting">
<ns0:netcdf location="http://gws-access.ceda.ac.uk/public/nemo/runs/ORCA0083-N01/means/1979/ORCA0083-N01_19791206d05T.nc" />
</ns0:aggregation>
</ns0:netcdf>
<ns0:netcdf>
<ns0:aggregation dimName="time_counter" name="zonal_velocity" type="joinExisting">
<ns0:netcdf location="http://gws-access.ceda.ac.uk/public/nemo/runs/ORCA0083-N01/means/1979/ORCA0083-N01_19791206d05U.nc" />
</ns0:aggregation>
</ns0:netcdf>
<ns0:netcdf>
<ns0:aggregation dimName="time_counter" name="meridian_velocity" type="joinExisting">
<ns0:netcdf location="http://gws-access.ceda.ac.uk/public/nemo/runs/ORCA0083-N01/means/1979/ORCA0083-N01_19791206d05V.nc" />
</ns0:aggregation>
</ns0:netcdf>
<ns0:netcdf>
<ns0:aggregation dimName="time_counter" name="sea_surface_height" type="joinExisting">
<ns0:netcdf location="http://gws-access.ceda.ac.uk/public/nemo/runs/ORCA0083-N01/means/1979/ORCA0083-N01_19791206d05T.nc" />
</ns0:aggregation>
</ns0:netcdf>
</ns0:aggregation>
</ns0:netcdf>
NAMELISTS_AND_FORTRAN_FILES/TIDES/namelist.bdy 0 → 100755
View file @ 868014cc
!-----------------------------------------------------------------------
! vertical coordinate
!-----------------------------------------------------------------------
ln_zco = .true. ! z-coordinate - full steps (T/F)
ln_zps = .false. ! z-coordinate - partial steps (T/F)
ln_sco = .false. ! s- or hybrid z-s-coordinate (T/F)
rn_hmin = -5 ! min depth of the ocean (>0) or
! min number of ocean level (<0)
!-----------------------------------------------------------------------
! s-coordinate or hybrid z-s-coordinate
!-----------------------------------------------------------------------
rn_sbot_min = 10. ! minimum depth of s-bottom surface (>0) (m)
rn_sbot_max = 7000. ! maximum depth of s-bottom surface
! (= ocean depth) (>0) (m)
ln_s_sigma = .false. ! hybrid s-sigma coordinates
rn_hc = 50.0 ! critical depth with s-sigma
!-----------------------------------------------------------------------
! grid information
!-----------------------------------------------------------------------
sn_src_hgr = './mesh_hgr_src.nc' ! parent /grid/
sn_src_zgr = './mesh_zgr_src.nc' ! parent
sn_dst_hgr = './domain_cfg.nc'
sn_dst_zgr = './inputs_dst.ncml' ! rename output variables
sn_src_msk = './mask_src.nc' ! parent
sn_bathy = './bathy_meter.nc'
!-----------------------------------------------------------------------
! I/O
!-----------------------------------------------------------------------
sn_src_dir = './inputs_src.ncml' ! src_files/'
sn_dst_dir = './'
sn_fn = 'INDIAN' ! prefix for output files
nn_fv = -1e20 ! set fill value for output files
nn_src_time_adj = 0 ! src time adjustment
sn_dst_metainfo = 'metadata info: ashbre'
!-----------------------------------------------------------------------
! unstructured open boundaries
!-----------------------------------------------------------------------
ln_coords_file = .true. ! =T : produce bdy coordinates files
cn_coords_file = 'coordinates.bdy.nc' ! name of bdy coordinates files (if ln_coords_file=.TRUE.)
ln_mask_file = .false. ! =T : read mask from file
cn_mask_file = './bdy_mask.nc' ! name of mask file (if ln_mask_file=.TRUE.)
ln_dyn2d = .true. ! boundary conditions for barotropic fields
ln_dyn3d = .false. ! boundary conditions for baroclinic velocities
ln_tra = .false. ! boundary conditions for T and S
ln_ice = .false. ! ice boundary condition
nn_rimwidth = 1 ! width of the relaxation zone
!-----------------------------------------------------------------------
! unstructured open boundaries tidal parameters
!-----------------------------------------------------------------------
ln_tide = .true. ! =T : produce bdy tidal conditions
clname(1) = 'M2'
clname(2) = 'S2'
clname(3) = 'K2'
clname(4) = 'Q1' ! name of constituent
clname(5) = 'O1'
clname(6) = 'P1'
clname(7) = 'K1'
clname(8) = 'N2'
ln_trans = .false.
sn_tide_h = './h_tpxo7.2.nc'
sn_tide_u = './u_tpxo7.2.nc'
!-----------------------------------------------------------------------
! Time information
!-----------------------------------------------------------------------
nn_year_000 = 1979 ! year start
nn_year_end = 1979 ! year end
nn_month_000 = 11 ! month start (default = 1 is years>1)
nn_month_end = 11 ! month end (default = 12 is years>1)
sn_dst_calendar = 'gregorian' ! output calendar format
nn_base_year = 1978 ! base year for time counter
sn_tide_grid = './grid_tpxo7.2.nc'
!-----------------------------------------------------------------------
! Additional parameters
!-----------------------------------------------------------------------
nn_wei = 1 ! smoothing filter weights
rn_r0 = 0.041666666 ! decorrelation distance use in gauss
! smoothing onto dst points. Need to
! make this a funct. of dlon
sn_history = 'bdy files produced by jelt from ORCA0083-N01'
! history for netcdf file
ln_nemo3p4 = .true. ! else presume v3.2 or v3.3
nn_alpha = 0 ! Euler rotation angle
nn_beta = 0 ! Euler rotation angle
nn_gamma = 0 ! Euler rotation angle
rn_mask_max_depth = 7000.0 ! Maximum depth to be ignored for the mask
rn_mask_shelfbreak_dist = 60 ! Distance from the shelf break
NAMELISTS_AND_FORTRAN_FILES/TIDES/run_script.pbs 0 → 100755
View file @ 868014cc
#!/bin/bash
# ---------------------------
#===============================================================
# CLUSTER BITS
#===============================================================
#PBS -N EA_R12
#PBS -l select=serial=true:ncpus=1
#PBS -l walltime=12:00:00
#PBS -A n01-ACCORD
#PBS -j oe
cd $PBS_O_WORKDIR
module unload anaconda
module load anaconda/2.2.0-python2
LD_LIBRARY_PATH=/opt/java/jdk1.8.0_51/jre/lib/amd64/server/:$LD_LIBRARY_PATH
export PYTHONPATH=~/.conda/envs/nrct_tide/lib/python2.7/site-packages/:$PYTHONPATH
source activate nrct_tide
#===============================================================
# LAUNCH JOB
#===============================================================
echo `date` : Launch Job
pynemo -s namelist.bdy
exit
NAMELISTS_AND_FORTRAN_FILES/f_files/bdyini.F90 0 → 100755
View file @ 868014cc
MODULE bdyini
!!======================================================================
!! *** MODULE bdyini ***
!! Unstructured open boundaries : initialisation
!!======================================================================
!! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code
!! - ! 2007-01 (D. Storkey) Update to use IOM module
!! - ! 2007-01 (D. Storkey) Tidal forcing
!! 3.0 ! 2008-04 (NEMO team) add in the reference version
!! 3.3 ! 2010-09 (E.O'Dea) updates for Shelf configurations
!! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions
!! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge
!! 3.4 ! 2012 (J. Chanut) straight open boundary case update
!! 3.5 ! 2012 (S. Mocavero, I. Epicoco) optimization of BDY communications
!! 3.7 ! 2016 (T. Lovato) Remove bdy macro, call here init for dta and tides
!!----------------------------------------------------------------------
!! bdy_init : Initialization of unstructured open boundaries
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE dom_oce ! ocean space and time domain
USE bdy_oce ! unstructured open boundary conditions
USE bdydta ! open boundary cond. setting (bdy_dta_init routine)
USE bdytides ! open boundary cond. setting (bdytide_init routine)
USE sbctide ! Tidal forcing or not
USE phycst , ONLY: rday
!
USE in_out_manager ! I/O units
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp ! for mpp_sum
USE iom ! I/O
USE wrk_nemo ! Memory Allocation
USE timing ! Timing
IMPLICIT NONE
PRIVATE
PUBLIC bdy_init ! routine called in nemo_init
INTEGER, PARAMETER :: jp_nseg = 100 !
INTEGER, PARAMETER :: nrimmax = 20 ! maximum rimwidth in structured
! open boundary data files
! Straight open boundary segment parameters:
INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs
INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge !
INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw !
INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn !
INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs !
!!----------------------------------------------------------------------
!! NEMO/OPA 3.7 , NEMO Consortium (2015)
!! $Id: bdyini.F90 7646 2017-02-06 09:25:03Z timgraham $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE bdy_init
!!----------------------------------------------------------------------
!! *** ROUTINE bdy_init ***
!!
!! ** Purpose : Initialization of the dynamics and tracer fields with
!! unstructured open boundaries.
!!
!! ** Method : Read initialization arrays (mask, indices) to identify
!! an unstructured open boundary
!!
!! ** Input : bdy_init.nc, input file for unstructured open boundaries
!!----------------------------------------------------------------------
NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file, &
& ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
& cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
& ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
& cn_ice_lim, nn_ice_lim_dta, &
& rn_ice_tem, rn_ice_sal, rn_ice_age, &
& ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
!
INTEGER :: ios ! Local integer output status for namelist read
!!----------------------------------------------------------------------
!
IF( nn_timing == 1 ) CALL timing_start('bdy_init')
! ------------------------
! Read namelist parameters
! ------------------------
REWIND( numnam_ref ) ! Namelist nambdy in reference namelist :Unstructured open boundaries
READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
!
REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist :Unstructured open boundaries
READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
IF(lwm) WRITE ( numond, nambdy )
! -----------------------------------------
! unstructured open boundaries use control
! -----------------------------------------
IF ( ln_bdy ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
IF(lwp) WRITE(numout,*) '~~~~~~~~'
!
! Open boundaries definition (arrays and masks)
CALL bdy_segs
!
! Open boundaries initialisation of external data arrays
CALL bdy_dta_init
!
! Open boundaries initialisation of tidal harmonic forcing
IF( ln_tide ) CALL bdytide_init
!
ELSE
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'bdy_init : open boundaries not used (ln_bdy = F)'
IF(lwp) WRITE(numout,*) '~~~~~~~~'
!
ENDIF
!
IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
!
END SUBROUTINE bdy_init
SUBROUTINE bdy_segs
!!----------------------------------------------------------------------
!! *** ROUTINE bdy_init ***
!!
!! ** Purpose : Definition of unstructured open boundaries.
!!
!! ** Method : Read initialization arrays (mask, indices) to identify
!! an unstructured open boundary
!!
!! ** Input : bdy_init.nc, input file for unstructured open boundaries
!!----------------------------------------------------------------------
! local variables
!-------------------
INTEGER :: ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
INTEGER :: icount, icountr, ibr_max, ilen1, ibm1 ! local integers
INTEGER :: iwe, ies, iso, ino, inum, id_dummy ! - -
INTEGER :: igrd_start, igrd_end, jpbdta ! - -
INTEGER :: jpbdtau, jpbdtas ! - -
INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - -
INTEGER :: i_offset, j_offset ! - -
INTEGER , POINTER :: nbi, nbj, nbr ! short cuts
REAL(wp), POINTER :: flagu, flagv ! - -
REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields
REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars
INTEGER, DIMENSION (2) :: kdimsz
INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays
INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta
INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points
CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid
INTEGER :: com_east, com_west, com_south, com_north ! Flags for boundaries sending
INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b ! Flags for boundaries receiving
INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4) ! Arrays for neighbours coordinates
REAL(wp), POINTER, DIMENSION(:,:) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat)
!!
CHARACTER(LEN=1) :: ctypebdy ! - -
INTEGER :: nbdyind, nbdybeg, nbdyend
!!
NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
INTEGER :: ios ! Local integer output status for namelist read
!!----------------------------------------------------------------------
!
IF( nn_timing == 1 ) CALL timing_start('bdy_segs')
!
cgrid = (/'t','u','v'/)
! -----------------------------------------
! Check and write out namelist parameters
! -----------------------------------------
! IF( jperio /= 0 ) CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,', &
! & ' and general open boundary condition are not compatible' )
IF( nb_bdy == 0 ) THEN
IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
ELSE
IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
ENDIF
DO ib_bdy = 1,nb_bdy
IF(lwp) WRITE(numout,*) ' '
IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'
IF( ln_coords_file(ib_bdy) ) THEN
IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
ELSE
IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
ENDIF
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: '
SELECT CASE( cn_dyn2d(ib_bdy) )
CASE( 'none' )
IF(lwp) WRITE(numout,*) ' no open boundary condition'
dta_bdy(ib_bdy)%ll_ssh = .false.
dta_bdy(ib_bdy)%ll_u2d = .false.
dta_bdy(ib_bdy)%ll_v2d = .false.
CASE( 'frs' )
IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
dta_bdy(ib_bdy)%ll_ssh = .false.
dta_bdy(ib_bdy)%ll_u2d = .true.
dta_bdy(ib_bdy)%ll_v2d = .true.
CASE( 'flather' )
IF(lwp) WRITE(numout,*) ' Flather radiation condition'
dta_bdy(ib_bdy)%ll_ssh = .true.
dta_bdy(ib_bdy)%ll_u2d = .true.
dta_bdy(ib_bdy)%ll_v2d = .true.
CASE( 'orlanski' )
IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
dta_bdy(ib_bdy)%ll_ssh = .false.
dta_bdy(ib_bdy)%ll_u2d = .true.
dta_bdy(ib_bdy)%ll_v2d = .true.
CASE( 'orlanski_npo' )
IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
dta_bdy(ib_bdy)%ll_ssh = .false.
dta_bdy(ib_bdy)%ll_u2d = .true.
dta_bdy(ib_bdy)%ll_v2d = .true.
CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
END SELECT
IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
SELECT CASE( nn_dyn2d_dta(ib_bdy) ) !
CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' tidal harmonic forcing taken from file'
CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files'
CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
END SELECT
IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.ln_tide)) THEN
CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )
ENDIF
ENDIF
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: '
SELECT CASE( cn_dyn3d(ib_bdy) )
CASE('none')
IF(lwp) WRITE(numout,*) ' no open boundary condition'
dta_bdy(ib_bdy)%ll_u3d = .false.
dta_bdy(ib_bdy)%ll_v3d = .false.
CASE('frs')
IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
dta_bdy(ib_bdy)%ll_u3d = .true.
dta_bdy(ib_bdy)%ll_v3d = .true.
CASE('specified')
IF(lwp) WRITE(numout,*) ' Specified value'
dta_bdy(ib_bdy)%ll_u3d = .true.
dta_bdy(ib_bdy)%ll_v3d = .true.
CASE('neumann')
IF(lwp) WRITE(numout,*) ' Neumann conditions'
dta_bdy(ib_bdy)%ll_u3d = .false.
dta_bdy(ib_bdy)%ll_v3d = .false.
CASE('zerograd')
IF(lwp) WRITE(numout,*) ' Zero gradient for baroclinic velocities'
dta_bdy(ib_bdy)%ll_u3d = .false.
dta_bdy(ib_bdy)%ll_v3d = .false.
CASE('zero')
IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)'
dta_bdy(ib_bdy)%ll_u3d = .false.
dta_bdy(ib_bdy)%ll_v3d = .false.
CASE('orlanski')
IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
dta_bdy(ib_bdy)%ll_u3d = .true.
dta_bdy(ib_bdy)%ll_v3d = .true.
CASE('orlanski_npo')
IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
dta_bdy(ib_bdy)%ll_u3d = .true.
dta_bdy(ib_bdy)%ll_v3d = .true.
CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
END SELECT
IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
SELECT CASE( nn_dyn3d_dta(ib_bdy) ) !
CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
END SELECT
ENDIF
IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
ln_dyn3d_dmp(ib_bdy)=.false.
ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
CALL ctl_stop( 'Use FRS OR relaxation' )
ELSE
IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone'
IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'
IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
dta_bdy(ib_bdy)%ll_u3d = .true.
dta_bdy(ib_bdy)%ll_v3d = .true.
ENDIF
ELSE
IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities'
ENDIF
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: '
SELECT CASE( cn_tra(ib_bdy) )
CASE('none')
IF(lwp) WRITE(numout,*) ' no open boundary condition'
dta_bdy(ib_bdy)%ll_tem = .false.
dta_bdy(ib_bdy)%ll_sal = .false.
CASE('frs')
IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
dta_bdy(ib_bdy)%ll_tem = .true.
dta_bdy(ib_bdy)%ll_sal = .true.
CASE('specified')
IF(lwp) WRITE(numout,*) ' Specified value'
dta_bdy(ib_bdy)%ll_tem = .true.
dta_bdy(ib_bdy)%ll_sal = .true.
CASE('neumann')
IF(lwp) WRITE(numout,*) ' Neumann conditions'
dta_bdy(ib_bdy)%ll_tem = .false.
dta_bdy(ib_bdy)%ll_sal = .false.
CASE('runoff')
IF(lwp) WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity'
dta_bdy(ib_bdy)%ll_tem = .false.
dta_bdy(ib_bdy)%ll_sal = .false.
CASE('orlanski')
IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
dta_bdy(ib_bdy)%ll_tem = .true.
dta_bdy(ib_bdy)%ll_sal = .true.
CASE('orlanski_npo')
IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
dta_bdy(ib_bdy)%ll_tem = .true.
dta_bdy(ib_bdy)%ll_sal = .true.
CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_tra' )
END SELECT
IF( cn_tra(ib_bdy) /= 'none' ) THEN
SELECT CASE( nn_tra_dta(ib_bdy) ) !
CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
END SELECT
ENDIF
IF ( ln_tra_dmp(ib_bdy) ) THEN
IF ( cn_tra(ib_bdy) == 'none' ) THEN
IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
ln_tra_dmp(ib_bdy)=.false.
ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
CALL ctl_stop( 'Use FRS OR relaxation' )
ELSE
IF(lwp) WRITE(numout,*) ' + T/S relaxation zone'
IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'
IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
dta_bdy(ib_bdy)%ll_tem = .true.
dta_bdy(ib_bdy)%ll_sal = .true.
ENDIF
ELSE
IF(lwp) WRITE(numout,*) ' NO T/S relaxation'
ENDIF
IF(lwp) WRITE(numout,*)
#if defined key_lim2
IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: '
SELECT CASE( cn_ice_lim(ib_bdy) )
CASE('none')
IF(lwp) WRITE(numout,*) ' no open boundary condition'
dta_bdy(ib_bdy)%ll_frld = .false.
dta_bdy(ib_bdy)%ll_hicif = .false.
dta_bdy(ib_bdy)%ll_hsnif = .false.
CASE('frs')
IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
dta_bdy(ib_bdy)%ll_frld = .true.
dta_bdy(ib_bdy)%ll_hicif = .true.
dta_bdy(ib_bdy)%ll_hsnif = .true.
CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice_lim' )
END SELECT
IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN
SELECT CASE( nn_ice_lim_dta(ib_bdy) ) !
CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' )
END SELECT
ENDIF
IF(lwp) WRITE(numout,*)
#elif defined key_lim3
IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: '
SELECT CASE( cn_ice_lim(ib_bdy) )
CASE('none')
IF(lwp) WRITE(numout,*) ' no open boundary condition'
dta_bdy(ib_bdy)%ll_a_i = .false.
dta_bdy(ib_bdy)%ll_ht_i = .false.
dta_bdy(ib_bdy)%ll_ht_s = .false.
CASE('frs')
IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
dta_bdy(ib_bdy)%ll_a_i = .true.
dta_bdy(ib_bdy)%ll_ht_i = .true.
dta_bdy(ib_bdy)%ll_ht_s = .true.
CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice_lim' )
END SELECT
IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN
SELECT CASE( nn_ice_lim_dta(ib_bdy) ) !
CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'
CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'
CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' )
END SELECT
ENDIF
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)
IF(lwp) WRITE(numout,*) ' sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)
IF(lwp) WRITE(numout,*) ' age of bdy sea-ice = ', rn_ice_age(ib_bdy)
#endif
IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy)
IF(lwp) WRITE(numout,*)
ENDDO
IF (nb_bdy .gt. 0) THEN
IF( ln_vol ) THEN ! check volume conservation (nn_volctl value)
IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
IF(lwp) WRITE(numout,*)
SELECT CASE ( nn_volctl )
CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant'
CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux'
CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' )
END SELECT
IF(lwp) WRITE(numout,*)
ELSE
IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
IF(lwp) WRITE(numout,*)
ENDIF
IF( nb_jpk_bdy > 0 ) THEN
IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'
ELSE
IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***'
ENDIF
ENDIF
! -------------------------------------------------
! Initialise indices arrays for open boundaries
! -------------------------------------------------
! Work out global dimensions of boundary data
! ---------------------------------------------
REWIND( numnam_cfg )
nblendta(:,:) = 0
nbdysege = 0
nbdysegw = 0
nbdysegn = 0
nbdysegs = 0
icount = 0 ! count user defined segments
! Dimensions below are used to allocate arrays to read external data
jpbdtas = 1 ! Maximum size of boundary data (structured case)
jpbdtau = 1 ! Maximum size of boundary data (unstructured case)
DO ib_bdy = 1, nb_bdy
IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
icount = icount + 1
! No REWIND here because may need to read more than one nambdy_index namelist.
! Read only namelist_cfg to avoid unseccessfull overwrite
!! REWIND( numnam_ref ) ! Namelist nambdy_index in reference namelist : Open boundaries indexes
!! READ ( numnam_ref, namrun, IOSTAT = ios, ERR = 903)
!!903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in reference namelist', lwp )
!! REWIND( numnam_cfg ) ! Namelist nambdy_index in configuration namelist : Open boundaries indexes
READ ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp )
IF(lwm) WRITE ( numond, nambdy_index )
SELECT CASE ( TRIM(ctypebdy) )
CASE( 'N' )
IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
nbdyind = jpjglo - 2 ! set boundary to whole side of model domain.
nbdybeg = 2
nbdyend = jpiglo - 1
ENDIF
nbdysegn = nbdysegn + 1
npckgn(nbdysegn) = ib_bdy ! Save bdy package number
jpjnob(nbdysegn) = nbdyind
jpindt(nbdysegn) = nbdybeg
jpinft(nbdysegn) = nbdyend
!
CASE( 'S' )
IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
nbdyind = 2 ! set boundary to whole side of model domain.
nbdybeg = 2
nbdyend = jpiglo - 1
ENDIF
nbdysegs = nbdysegs + 1
npckgs(nbdysegs) = ib_bdy ! Save bdy package number
jpjsob(nbdysegs) = nbdyind
jpisdt(nbdysegs) = nbdybeg
jpisft(nbdysegs) = nbdyend
!
CASE( 'E' )
IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
nbdyind = jpiglo - 2 ! set boundary to whole side of model domain.
nbdybeg = 2
nbdyend = jpjglo - 1
ENDIF
nbdysege = nbdysege + 1
npckge(nbdysege) = ib_bdy ! Save bdy package number
jpieob(nbdysege) = nbdyind
jpjedt(nbdysege) = nbdybeg
jpjeft(nbdysege) = nbdyend
!
CASE( 'W' )
IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
nbdyind = 2 ! set boundary to whole side of model domain.
nbdybeg = 2
nbdyend = jpjglo - 1
ENDIF
nbdysegw = nbdysegw + 1
npckgw(nbdysegw) = ib_bdy ! Save bdy package number
jpiwob(nbdysegw) = nbdyind
jpjwdt(nbdysegw) = nbdybeg
jpjwft(nbdysegw) = nbdyend
!
CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
END SELECT
! For simplicity we assume that in case of straight bdy, arrays have the same length
! (even if it is true that last tangential velocity points
! are useless). This simplifies a little bit boundary data format (and agrees with format
! used so far in obc package)
nblendta(1:jpbgrd,ib_bdy) = (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy)
jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1))
IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) &
& CALL ctl_stop( 'rimwidth must be lower than nrimmax' )
ELSE ! Read size of arrays in boundary coordinates file.
CALL iom_open( cn_coords_file(ib_bdy), inum )
DO igrd = 1, jpbgrd
id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )
!clem nblendta(igrd,ib_bdy) = kdimsz(1)
!clem jpbdtau = MAX(jpbdtau, kdimsz(1))
nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))
END DO
CALL iom_close( inum )
!
ENDIF
!
END DO ! ib_bdy
IF (nb_bdy>0) THEN
jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
! Allocate arrays
!---------------
ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), &
& nbrdta(jpbdta, jpbgrd, nb_bdy) )
IF( nb_jpk_bdy>0 ) THEN
ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) )
ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) )
ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) )
ELSE
ALLOCATE( dta_global(jpbdtau, 1, jpk) )
ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO
ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO
ENDIF
IF ( icount>0 ) THEN
IF( nb_jpk_bdy>0 ) THEN
ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) )
ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) )
ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) )
ELSE
ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO
ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO
ENDIF
ENDIF
!
ENDIF
! Now look for crossings in user (namelist) defined open boundary segments:
!--------------------------------------------------------------------------
IF( icount>0 ) CALL bdy_ctl_seg
! Calculate global boundary index arrays or read in from file
!------------------------------------------------------------
! 1. Read global index arrays from boundary coordinates file.
DO ib_bdy = 1, nb_bdy
!
IF( ln_coords_file(ib_bdy) ) THEN
!
CALL iom_open( cn_coords_file(ib_bdy), inum )
DO igrd = 1, jpbgrd
CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
DO ii = 1,nblendta(igrd,ib_bdy)
nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
END DO
CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
DO ii = 1,nblendta(igrd,ib_bdy)
nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
END DO
CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
DO ii = 1,nblendta(igrd,ib_bdy)
nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
END DO
!
ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
IF (ibr_max < nn_rimwidth(ib_bdy)) &
CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
END DO
CALL iom_close( inum )
!
ENDIF
!
END DO
! 2. Now fill indices corresponding to straight open boundary arrays:
! East
!-----
DO iseg = 1, nbdysege
ib_bdy = npckge(iseg)
!
! ------------ T points -------------
igrd=1
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
DO ij = jpjedt(iseg), jpjeft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
nbjdta(icount, igrd, ib_bdy) = ij
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
ENDDO
!
! ------------ U points -------------
igrd=2
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
DO ij = jpjedt(iseg), jpjeft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
nbjdta(icount, igrd, ib_bdy) = ij
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
ENDDO
!
! ------------ V points -------------
igrd=3
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
! DO ij = jpjedt(iseg), jpjeft(iseg) - 1
DO ij = jpjedt(iseg), jpjeft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
nbjdta(icount, igrd, ib_bdy) = ij
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
ENDDO
ENDDO
!
! West
!-----
DO iseg = 1, nbdysegw
ib_bdy = npckgw(iseg)
!
! ------------ T points -------------
igrd=1
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
DO ij = jpjwdt(iseg), jpjwft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
nbjdta(icount, igrd, ib_bdy) = ij
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
ENDDO
!
! ------------ U points -------------
igrd=2
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
DO ij = jpjwdt(iseg), jpjwft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
nbjdta(icount, igrd, ib_bdy) = ij
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
ENDDO
!
! ------------ V points -------------
igrd=3
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
DO ij = jpjwdt(iseg), jpjwft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
nbjdta(icount, igrd, ib_bdy) = ij
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
ENDDO
ENDDO
!
! North
!-----
DO iseg = 1, nbdysegn
ib_bdy = npckgn(iseg)
!
! ------------ T points -------------
igrd=1
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
DO ii = jpindt(iseg), jpinft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = ii
nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
ENDDO
!
! ------------ U points -------------
igrd=2
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
! DO ii = jpindt(iseg), jpinft(iseg) - 1
DO ii = jpindt(iseg), jpinft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = ii
nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
ENDDO
!
! ------------ V points -------------
igrd=3
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
DO ii = jpindt(iseg), jpinft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = ii
nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
ENDDO
ENDDO
!
! South
!-----
DO iseg = 1, nbdysegs
ib_bdy = npckgs(iseg)
!
! ------------ T points -------------
igrd=1
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
DO ii = jpisdt(iseg), jpisft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = ii
nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
ENDDO
!
! ------------ U points -------------
igrd=2
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
! DO ii = jpisdt(iseg), jpisft(iseg) - 1
DO ii = jpisdt(iseg), jpisft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = ii
nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
ENDDO
!
! ------------ V points -------------
igrd=3
icount=0
DO ir = 1, nn_rimwidth(ib_bdy)
DO ii = jpisdt(iseg), jpisft(iseg)
icount = icount + 1
nbidta(icount, igrd, ib_bdy) = ii
nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
nbrdta(icount, igrd, ib_bdy) = ir
ENDDO
ENDDO
ENDDO
! Deal with duplicated points
!-----------------------------
! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
! if their distance to the bdy is greater than the other
! If their distance are the same, just keep only one to avoid updating a point twice
DO igrd = 1, jpbgrd
DO ib_bdy1 = 1, nb_bdy
DO ib_bdy2 = 1, nb_bdy
IF (ib_bdy1/=ib_bdy2) THEN
DO ib1 = 1, nblendta(igrd,ib_bdy1)
DO ib2 = 1, nblendta(igrd,ib_bdy2)
IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
& (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
! & nbidta(ib1, igrd, ib_bdy1), &
! & nbjdta(ib2, igrd, ib_bdy2)
! keep only points with the lowest distance to boundary:
IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
! Arbitrary choice if distances are the same:
ELSE
nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
ENDIF
END IF
END DO
END DO
ENDIF
END DO
END DO
END DO
! Work out dimensions of boundary data on each processor
! ------------------------------------------------------
! Rather assume that boundary data indices are given on global domain
! TO BE DISCUSSED ?
! iw = mig(1) + 1 ! if monotasking and no zoom, iw=2
! ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
! is = mjg(1) + 1 ! if monotasking and no zoom, is=2
! in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
iwe = mig(1) - 1 + 2 ! if monotasking and no zoom, iw=2
ies = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
iso = mjg(1) - 1 + 2 ! if monotasking and no zoom, is=2
ino = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
ALLOCATE( nbondi_bdy(nb_bdy))
ALLOCATE( nbondj_bdy(nb_bdy))
nbondi_bdy(:)=2
nbondj_bdy(:)=2
ALLOCATE( nbondi_bdy_b(nb_bdy))
ALLOCATE( nbondj_bdy_b(nb_bdy))
nbondi_bdy_b(:)=2
nbondj_bdy_b(:)=2
! Work out dimensions of boundary data on each neighbour process
IF(nbondi == 0) THEN
iw_b(1) = 1 + nimppt(nowe+1)
ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3
is_b(1) = 1 + njmppt(nowe+1)
in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3
iw_b(2) = 1 + nimppt(noea+1)
ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3
is_b(2) = 1 + njmppt(noea+1)
in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3
ELSEIF(nbondi == 1) THEN
iw_b(1) = 1 + nimppt(nowe+1)
ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3
is_b(1) = 1 + njmppt(nowe+1)
in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3
ELSEIF(nbondi == -1) THEN
iw_b(2) = 1 + nimppt(noea+1)
ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3
is_b(2) = 1 + njmppt(noea+1)
in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3
ENDIF
IF(nbondj == 0) THEN
iw_b(3) = 1 + nimppt(noso+1)
ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3
is_b(3) = 1 + njmppt(noso+1)
in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3
iw_b(4) = 1 + nimppt(nono+1)
ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3
is_b(4) = 1 + njmppt(nono+1)
in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3
ELSEIF(nbondj == 1) THEN
iw_b(3) = 1 + nimppt(noso+1)
ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3
is_b(3) = 1 + njmppt(noso+1)
in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3
ELSEIF(nbondj == -1) THEN
iw_b(4) = 1 + nimppt(nono+1)
ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3
is_b(4) = 1 + njmppt(nono+1)
in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3
ENDIF
DO ib_bdy = 1, nb_bdy
DO igrd = 1, jpbgrd
icount = 0
icountr = 0
idx_bdy(ib_bdy)%nblen(igrd) = 0
idx_bdy(ib_bdy)%nblenrim(igrd) = 0
DO ib = 1, nblendta(igrd,ib_bdy)
! check that data is in correct order in file
ibm1 = MAX(1,ib-1)
IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc...
IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
& ' in order of distance from edge nbr A utility for re-ordering ', &
& ' boundary coordinates and data files exists in the TOOLS/OBC directory')
ENDIF
ENDIF
! check if point is in local domain
IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
& nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino ) THEN
!
icount = icount + 1
!
IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr+1
ENDIF
ENDDO
idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc
ENDDO ! igrd
! Allocate index arrays for this boundary set
!--------------------------------------------
ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
ALLOCATE( idx_bdy(ib_bdy)%nbi (ilen1,jpbgrd) )
ALLOCATE( idx_bdy(ib_bdy)%nbj (ilen1,jpbgrd) )
ALLOCATE( idx_bdy(ib_bdy)%nbr (ilen1,jpbgrd) )
ALLOCATE( idx_bdy(ib_bdy)%nbd (ilen1,jpbgrd) )
ALLOCATE( idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) )
ALLOCATE( idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) )
ALLOCATE( idx_bdy(ib_bdy)%nbw (ilen1,jpbgrd) )
ALLOCATE( idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) )
ALLOCATE( idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
! Dispatch mapping indices and discrete distances on each processor
! -----------------------------------------------------------------
com_east = 0
com_west = 0
com_south = 0
com_north = 0
com_east_b = 0
com_west_b = 0
com_south_b = 0
com_north_b = 0
DO igrd = 1, jpbgrd
icount = 0
! Loop on rimwidth to ensure outermost points come first in the local arrays.
DO ir=1, nn_rimwidth(ib_bdy)
DO ib = 1, nblendta(igrd,ib_bdy)
! check if point is in local domain and equals ir
IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
& nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
!
icount = icount + 1
! Rather assume that boundary data indices are given on global domain
! TO BE DISCUSSED ?
! idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
! idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
! check if point has to be sent
ii = idx_bdy(ib_bdy)%nbi(icount,igrd)
ij = idx_bdy(ib_bdy)%nbj(icount,igrd)
if((com_east .ne. 1) .and. (ii == (nlci-1)) .and. (nbondi .le. 0)) then
com_east = 1
elseif((com_west .ne. 1) .and. (ii == 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
com_west = 1
endif
if((com_south .ne. 1) .and. (ij == 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
com_south = 1
elseif((com_north .ne. 1) .and. (ij == (nlcj-1)) .and. (nbondj .le. 0)) then
com_north = 1
endif
idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy)
idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
ENDIF
! check if point has to be received from a neighbour
IF(nbondi == 0) THEN
IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. &
& nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
com_south = 1
elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
com_north = 1
endif
com_west_b = 1
endif
ENDIF
IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. &
& nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
if((com_east_b .ne. 1) .and. (ii == 2)) then
ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
com_south = 1
elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
com_north = 1
endif
com_east_b = 1
endif
ENDIF
ELSEIF(nbondi == 1) THEN
IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. &
& nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
com_south = 1
elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
com_north = 1
endif
com_west_b = 1
endif
ENDIF
ELSEIF(nbondi == -1) THEN
IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. &
& nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
if((com_east_b .ne. 1) .and. (ii == 2)) then
ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
com_south = 1
elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
com_north = 1
endif
com_east_b = 1
endif
ENDIF
ENDIF
IF(nbondj == 0) THEN
IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 &
& .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
& nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
com_north_b = 1
ENDIF
IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 &
&.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
& nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
com_south_b = 1
ENDIF
IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. &
& nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then
com_south_b = 1
endif
ENDIF
IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. &
& nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
if((com_north_b .ne. 1) .and. (ij == 2)) then
com_north_b = 1
endif
ENDIF
ELSEIF(nbondj == 1) THEN
IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. &
& nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
& nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
com_south_b = 1
ENDIF
IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. &
& nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then
com_south_b = 1
endif
ENDIF
ELSEIF(nbondj == -1) THEN
IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 &
& .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
& nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
com_north_b = 1
ENDIF
IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. &
& nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. &
& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
if((com_north_b .ne. 1) .and. (ij == 2)) then
com_north_b = 1
endif
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
! definition of the i- and j- direction local boundaries arrays used for sending the boundaries
IF( (com_east == 1) .and. (com_west == 1) ) THEN ; nbondi_bdy(ib_bdy) = 0
ELSEIF( (com_east == 1) .and. (com_west == 0) ) THEN ; nbondi_bdy(ib_bdy) = -1
ELSEIF( (com_east == 0) .and. (com_west == 1) ) THEN ; nbondi_bdy(ib_bdy) = 1
ENDIF
IF( (com_north == 1) .and. (com_south == 1) ) THEN ; nbondj_bdy(ib_bdy) = 0
ELSEIF( (com_north == 1) .and. (com_south == 0) ) THEN ; nbondj_bdy(ib_bdy) = -1
ELSEIF( (com_north == 0) .and. (com_south == 1) ) THEN ; nbondj_bdy(ib_bdy) = 1
ENDIF
! definition of the i- and j- direction local boundaries arrays used for receiving the boundaries
IF( (com_east_b == 1) .and. (com_west_b == 1) ) THEN ; nbondi_bdy_b(ib_bdy) = 0
ELSEIF( (com_east_b == 1) .and. (com_west_b == 0) ) THEN ; nbondi_bdy_b(ib_bdy) = -1
ELSEIF( (com_east_b == 0) .and. (com_west_b == 1) ) THEN ; nbondi_bdy_b(ib_bdy) = 1
ENDIF
IF( (com_north_b == 1) .and. (com_south_b == 1) ) THEN ; nbondj_bdy_b(ib_bdy) = 0
ELSEIF( (com_north_b == 1) .and. (com_south_b == 0) ) THEN ; nbondj_bdy_b(ib_bdy) = -1
ELSEIF( (com_north_b == 0) .and. (com_south_b == 1) ) THEN ; nbondj_bdy_b(ib_bdy) = 1
ENDIF
! Compute rim weights for FRS scheme
! ----------------------------------
DO igrd = 1, jpbgrd
DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( nbr - 1 ) *0.5 ) ! tanh formulation
! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)) ! linear
END DO
END DO
! Compute damping coefficients
! ----------------------------
DO igrd = 1, jpbgrd
DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &
& *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &
& *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
END DO
END DO
ENDDO
! ------------------------------------------------------
! Initialise masks and find normal/tangential directions
! ------------------------------------------------------
! Read global 2D mask at T-points: bdytmask
! -----------------------------------------
! bdytmask = 1 on the computational domain AND on open boundaries
! = 0 elsewhere
bdytmask(:,:) = ssmask(:,:)
! we need to derive mask on U and V grid from mask on T grid here.
bdyumask(:,:) = 0._wp
bdyvmask(:,:) = 0._wp
DO ij = 1, jpjm1
DO ii = 1, jpim1
bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij )
bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii ,ij+1)
END DO
END DO
CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond.
! bdy masks are now set to zero on boundary points:
!
igrd = 1
DO ib_bdy = 1, nb_bdy
DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
END DO
END DO
!
igrd = 2
DO ib_bdy = 1, nb_bdy
DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
END DO
END DO
!
igrd = 3
DO ib_bdy = 1, nb_bdy
DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
ENDDO
ENDDO
! For the flagu/flagv calculation below we require a version of fmask without
! the land boundary condition (shlat) included:
CALL wrk_alloc(jpi,jpj, zfmask )
DO ij = 2, jpjm1
DO ii = 2, jpim1
zfmask(ii,ij) = tmask(ii,ij ,1) * tmask(ii+1,ij ,1) &
& * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1)
END DO
END DO
! Lateral boundary conditions
CALL lbc_lnk( zfmask , 'F', 1. )
CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
idx_bdy(ib_bdy)%flagu(:,:) = 0._wp
idx_bdy(ib_bdy)%flagv(:,:) = 0._wp
icount = 0
! Calculate relationship of U direction to the local orientation of the boundary
! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
! flagu = 0 : u is tangential
! flagu = 1 : u is normal to the boundary and is direction is inward
DO igrd = 1,jpbgrd
SELECT CASE( igrd )
CASE( 1 ) ; pmask => umask (:,:,1) ; i_offset = 0
CASE( 2 ) ; pmask => bdytmask(:,:) ; i_offset = 1
CASE( 3 ) ; pmask => zfmask (:,:) ; i_offset = 0
END SELECT
icount = 0
DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
zefl = pmask(nbi+i_offset-1,nbj)
zwfl = pmask(nbi+i_offset,nbj)
! This error check only works if you are using the bdyXmask arrays
IF( i_offset == 1 .and. zefl + zwfl == 2 ) THEN
icount = icount + 1
IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
ELSE
idx_bdy(ib_bdy)%flagu(ib,igrd) = -zefl + zwfl
ENDIF
END DO
IF( icount /= 0 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', &
' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
IF(lwp) WRITE(numout,*) ' ========== '
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ENDIF
END DO
! Calculate relationship of V direction to the local orientation of the boundary
! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
! flagv = 0 : v is tangential
! flagv = 1 : v is normal to the boundary and is direction is inward
DO igrd = 1, jpbgrd
SELECT CASE( igrd )
CASE( 1 ) ; pmask => vmask (:,:,1) ; j_offset = 0
CASE( 2 ) ; pmask => zfmask(:,:) ; j_offset = 0
CASE( 3 ) ; pmask => bdytmask ; j_offset = 1
END SELECT
icount = 0
DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
znfl = pmask(nbi,nbj+j_offset-1)
zsfl = pmask(nbi,nbj+j_offset )
! This error check only works if you are using the bdyXmask arrays
IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN
IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
icount = icount + 1
ELSE
idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl
END IF
END DO
IF( icount /= 0 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', &
' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
IF(lwp) WRITE(numout,*) ' ========== '
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ENDIF
END DO
!
END DO
! Compute total lateral surface for volume correction:
! ----------------------------------------------------
! JC: this must be done at each time step with non-linear free surface
bdysurftot = 0._wp
IF( ln_vol ) THEN
igrd = 2 ! Lateral surface at U-points
DO ib_bdy = 1, nb_bdy
DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
flagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
bdysurftot = bdysurftot + hu_n (nbi , nbj) &
& * e2u (nbi , nbj) * ABS( flagu ) &
& * tmask_i(nbi , nbj) &
& * tmask_i(nbi+1, nbj)
END DO
END DO
igrd=3 ! Add lateral surface at V-points
DO ib_bdy = 1, nb_bdy
DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
flagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
bdysurftot = bdysurftot + hv_n (nbi, nbj ) &
& * e1v (nbi, nbj ) * ABS( flagv ) &
& * tmask_i(nbi, nbj ) &
& * tmask_i(nbi, nbj+1)
END DO
END DO
!
IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain
END IF
!
! Tidy up
!--------
IF( nb_bdy>0 ) DEALLOCATE( nbidta, nbjdta, nbrdta )
!
CALL wrk_dealloc(jpi,jpj, zfmask )
!
IF( nn_timing == 1 ) CALL timing_stop('bdy_segs')
!
END SUBROUTINE bdy_segs
SUBROUTINE bdy_ctl_seg
!!----------------------------------------------------------------------
!! *** ROUTINE bdy_ctl_seg ***
!!
!! ** Purpose : Check straight open boundary segments location
!!
!! ** Method : - Look for open boundary corners
!! - Check that segments start or end on land
!!----------------------------------------------------------------------
INTEGER :: ib, ib1, ib2, ji ,jj, itest
INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns
REAL(wp), DIMENSION(2) :: ztestmask
!!----------------------------------------------------------------------
!
IF (lwp) WRITE(numout,*) ' '
IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
!
IF(lwp) WRITE(numout,*) 'Number of east segments : ', nbdysege
IF(lwp) WRITE(numout,*) 'Number of west segments : ', nbdysegw
IF(lwp) WRITE(numout,*) 'Number of north segments : ', nbdysegn
IF(lwp) WRITE(numout,*) 'Number of south segments : ', nbdysegs
! 1. Check bounds
!----------------
DO ib = 1, nbdysegn
IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
IF ((jpjnob(ib).ge.jpjglo-1).or.&
&(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
IF (jpindt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
IF (jpinft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' )
END DO
!
DO ib = 1, nbdysegs
IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
IF ((jpjsob(ib).ge.jpjglo-1).or.&
&(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
IF (jpisdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
IF (jpisft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' )
END DO
!
DO ib = 1, nbdysege
IF (lwp) WRITE(numout,*) '**check east seg bounds pckg: ', npckge(ib)
IF ((jpieob(ib).ge.jpiglo-1).or.&
&(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
IF (jpjedt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
IF (jpjeft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' )
END DO
!
DO ib = 1, nbdysegw
IF (lwp) WRITE(numout,*) '**check west seg bounds pckg: ', npckgw(ib)
IF ((jpiwob(ib).ge.jpiglo-1).or.&
&(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
IF (jpjwdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
IF (jpjwft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' )
ENDDO
!
!
! 2. Look for segment crossings
!------------------------------
IF (lwp) WRITE(numout,*) '**Look for segments corners :'
!
itest = 0 ! corner number
!
! flag to detect if start or end of open boundary belongs to a corner
! if not (=0), it must be on land.
! if a corner is detected, save bdy package number for further tests
icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
! South/West crossings
IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
DO ib1 = 1, nbdysegw
DO ib2 = 1, nbdysegs
IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
& ( jpisft(ib2)>=jpiwob(ib1)).AND. &
& ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
& ( jpjwft(ib1)>=jpjsob(ib2))) THEN
IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN
! We have a possible South-West corner
! WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
icornw(ib1,1) = npckgs(ib2)
icorns(ib2,1) = npckgw(ib1)
ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
& jpisft(ib2), jpjwft(ib1)
IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
IF(lwp) WRITE(numout,*) ' Crossing problem with West segment: ',npckgw(ib1), &
& ' and South segment: ',npckgs(ib2)
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ELSE
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices'
IF(lwp) WRITE(numout,*) ' ========== Crossing problem with West segment: ',npckgw(ib1) , &
& ' and South segment: ',npckgs(ib2)
IF(lwp) WRITE(numout,*)
nstop = nstop+1
END IF
END IF
END DO
END DO
END IF
!
! South/East crossings
IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
DO ib1 = 1, nbdysege
DO ib2 = 1, nbdysegs
IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
& ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
& ( jpjedt(ib1)<=jpjsob(ib2) ).AND. &
& ( jpjeft(ib1)>=jpjsob(ib2) )) THEN
IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
! We have a possible South-East corner
! WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
icorne(ib1,1) = npckgs(ib2)
icorns(ib2,2) = npckge(ib1)
ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
& jpisdt(ib2), jpjeft(ib1)
IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
IF(lwp) WRITE(numout,*) ' Crossing problem with East segment: ',npckge(ib1), &
& ' and South segment: ',npckgs(ib2)
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ELSE
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices'
IF(lwp) WRITE(numout,*) ' ========== Crossing problem with East segment: ',npckge(ib1), &
& ' and South segment: ',npckgs(ib2)
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
END IF
END IF
END DO
END DO
END IF
!
! North/West crossings
IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
DO ib1 = 1, nbdysegw
DO ib2 = 1, nbdysegn
IF (( jpindt(ib2)<=jpiwob(ib1) ).AND. &
& ( jpinft(ib2)>=jpiwob(ib1) ).AND. &
& ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
& ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
! We have a possible North-West corner
! WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
icornw(ib1,2) = npckgn(ib2)
icornn(ib2,1) = npckgw(ib1)
ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
& jpinft(ib2), jpjwdt(ib1)
IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
IF(lwp) WRITE(numout,*) ' Crossing problem with West segment: ',npckgw(ib1), &
& ' and North segment: ',npckgn(ib2)
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ELSE
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices'
IF(lwp) WRITE(numout,*) ' ========== Crossing problem with West segment: ',npckgw(ib1), &
& ' and North segment: ',npckgn(ib2)
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
END IF
END IF
END DO
END DO
END IF
!
! North/East crossings
IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
DO ib1 = 1, nbdysege
DO ib2 = 1, nbdysegn
IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
& ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
& ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
& ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
! We have a possible North-East corner
! WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
icorne(ib1,2) = npckgn(ib2)
icornn(ib2,2) = npckge(ib1)
ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
& jpindt(ib2), jpjedt(ib1)
IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
IF(lwp) WRITE(numout,*) ' Crossing problem with East segment: ',npckge(ib1), &
& ' and North segment: ',npckgn(ib2)
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ELSE
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices'
IF(lwp) WRITE(numout,*) ' ========== Crossing problem with East segment: ',npckge(ib1), &
& ' and North segment: ',npckgn(ib2)
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
END IF
END IF
END DO
END DO
END IF
!
! 3. Check if segment extremities are on land
!--------------------------------------------
!
! West segments
DO ib = 1, nbdysegw
! get mask at boundary extremities:
ztestmask(1:2)=0.
DO ji = 1, jpi
DO jj = 1, jpj
IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
& ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
& ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)
END DO
END DO
IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
IF (ztestmask(1)==1) THEN
IF (icornw(ib,1)==0) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ELSE
! This is a corner
IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
itest=itest+1
ENDIF
ENDIF
IF (ztestmask(2)==1) THEN
IF (icornw(ib,2)==0) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ELSE
! This is a corner
IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
itest=itest+1
ENDIF
ENDIF
END DO
!
! East segments
DO ib = 1, nbdysege
! get mask at boundary extremities:
ztestmask(1:2)=0.
DO ji = 1, jpi
DO jj = 1, jpj
IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
& ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
& ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)
END DO
END DO
IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
IF (ztestmask(1)==1) THEN
IF (icorne(ib,1)==0) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ELSE
! This is a corner
IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
itest=itest+1
ENDIF
ENDIF
IF (ztestmask(2)==1) THEN
IF (icorne(ib,2)==0) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ELSE
! This is a corner
IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
itest=itest+1
ENDIF
ENDIF
END DO
!
! South segments
DO ib = 1, nbdysegs
! get mask at boundary extremities:
ztestmask(1:2)=0.
DO ji = 1, jpi
DO jj = 1, jpj
IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
& ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
& ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)
END DO
END DO
IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ENDIF
IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ENDIF
END DO
!
! North segments
DO ib = 1, nbdysegn
! get mask at boundary extremities:
ztestmask(1:2)=0.
DO ji = 1, jpi
DO jj = 1, jpj
IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
& ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
& ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)
END DO
END DO
IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
IF(lwp) WRITE(numout,*) ' ========== does not start on land'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ENDIF
IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
IF(lwp) WRITE(numout,*) ' ========== does not end on land'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ENDIF
END DO
!
IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
!
! Other tests TBD:
! segments completly on land
! optimized open boundary array length according to landmask
! Nudging layers that overlap with interior domain
!
END SUBROUTINE bdy_ctl_seg
SUBROUTINE bdy_ctl_corn( ib1, ib2 )
!!----------------------------------------------------------------------
!! *** ROUTINE bdy_ctl_corn ***
!!
!! ** Purpose : Check numerical schemes consistency between
!! segments having a common corner
!!
!! ** Method :
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: ib1, ib2
INTEGER :: itest
!!----------------------------------------------------------------------
itest = 0
IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) ) itest = itest + 1
IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) ) itest = itest + 1
IF( cn_tra (ib1) /= cn_tra (ib2) ) itest = itest + 1
!
IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) ) itest = itest + 1
IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) ) itest = itest + 1
IF( nn_tra_dta (ib1) /= nn_tra_dta (ib2) ) itest = itest + 1
!
IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) ) itest = itest + 1
!
IF( itest>0 ) THEN
IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
IF(lwp) WRITE(numout,*) ' ========== have different open bdy schemes'
IF(lwp) WRITE(numout,*)
nstop = nstop + 1
ENDIF
!
END SUBROUTINE bdy_ctl_corn
!!=================================================================================
END MODULE bdyini
NAMELISTS_AND_FORTRAN_FILES/f_files/diaharm.F90 0 → 100755
View file @ 868014cc
MODULE diaharm
!!======================================================================
!! *** MODULE diaharm ***
!! Harmonic analysis of tidal constituents
!!======================================================================
!! History : 3.1 ! 2007 (O. Le Galloudec, J. Chanut) Original code
!!----------------------------------------------------------------------
#if defined key_diaharm
!!----------------------------------------------------------------------
!! 'key_diaharm'
!!
!! NB: 2017-12 : add 3D harmonic analysis of velocities
!! integration of Maria Luneva's development
!! 'key_3Ddiaharm'
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE dom_oce ! ocean space and time domain
USE phycst
USE daymod
USE tide_mod
USE sbctide ! Tidal forcing or not
!
# if defined key_3Ddiaharm
USE zdf_oce
#endif
!
USE in_out_manager ! I/O units
USE iom ! I/0 library
USE ioipsl ! NetCDF IPSL library
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE timing ! preformance summary
USE wrk_nemo ! working arrays
IMPLICIT NONE
PRIVATE
LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .TRUE.
INTEGER, PARAMETER :: jpincomax = 2.*jpmax_harmo
INTEGER, PARAMETER :: jpdimsparse = jpincomax*300*24
! !!** namelist variables **
INTEGER :: nit000_han ! First time step used for harmonic analysis
INTEGER :: nitend_han ! Last time step used for harmonic analysis
INTEGER :: nstep_han ! Time step frequency for harmonic analysis
INTEGER :: nb_ana ! Number of harmonics to analyse
INTEGER , ALLOCATABLE, DIMENSION(:) :: name
REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut , vt , ft
# if defined key_3Ddiaharm
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:) :: ana_temp
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: out_eta , out_u, out_v , out_w , out_dzi
# else
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta , out_u, out_v
# endif
INTEGER :: ninco, nsparse
INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse
INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1
REAL(wp), DIMENSION(jpdimsparse) :: valuesparse
REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7
REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot
CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
PUBLIC dia_harm ! routine called by step.F90
!!----------------------------------------------------------------------
!! NEMO/OPA 3.5 , NEMO Consortium (2013)
!! $Id: diaharm.F90 5585 2015-07-10 14:19:11Z jchanut $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dia_harm_init
!!----------------------------------------------------------------------
!! *** ROUTINE dia_harm_init ***
!!
!! ** Purpose : Initialization of tidal harmonic analysis
!!
!! ** Method : Initialize frequency array and nodal factor for nit000_han
!!
!!--------------------------------------------------------------------
INTEGER :: jh, nhan, jl
INTEGER :: ios ! Local integer output status for namelist read
NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
!!----------------------------------------------------------------------
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
# if defined key_3Ddiaharm
WRITE(numout,*) ' - 3D harmonic analysis of currents actovated (key_3Ddiaharm)'
#endif
WRITE(numout,*) '~~~~~~~ '
ENDIF
!
IF( .NOT. ln_tide ) CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
!
CALL tide_init_Wave
!
REWIND( numnam_ref ) ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
READ ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in reference namelist', lwp )
REWIND( numnam_cfg ) ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis
READ ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist', lwp )
IF(lwm) WRITE ( numond, nam_diaharm )
!
IF(lwp) THEN
WRITE(numout,*) 'First time step used for analysis: nit000_han= ', nit000_han
WRITE(numout,*) 'Last time step used for analysis: nitend_han= ', nitend_han
WRITE(numout,*) 'Time step frequency for harmonic analysis: nstep_han= ', nstep_han
ENDIF
! Basic checks on harmonic analysis time window:
! ----------------------------------------------
IF( nit000 > nit000_han ) CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000', &
& ' restart capability not implemented' )
IF( nitend < nitend_han ) CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend', &
& 'restart capability not implemented' )
IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 ) &
& CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
nb_ana = 0
DO jh=1,jpmax_harmo
DO jl=1,jpmax_harmo
IF(TRIM(tname(jh)) == Wave(jl)%cname_tide) THEN
nb_ana=nb_ana+1
ENDIF
END DO
END DO
!
IF(lwp) THEN
WRITE(numout,*) ' Namelist nam_diaharm'
WRITE(numout,*) ' nb_ana = ', nb_ana
CALL flush(numout)
ENDIF
!
IF (nb_ana > jpmax_harmo) THEN
IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nb_ana must be lower than jpmax_harmo, stop'
IF(lwp) WRITE(numout,*) ' jpmax_harmo= ', jpmax_harmo
nstop = nstop + 1
ENDIF
ALLOCATE(name (nb_ana))
DO jh=1,nb_ana
DO jl=1,jpmax_harmo
IF (TRIM(tname(jh)) .eq. Wave(jl)%cname_tide) THEN
name(jh) = jl
EXIT
END IF
END DO
END DO
! Initialize frequency array:
! ---------------------------
ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
DO jh = 1, nb_ana
IF(lwp) WRITE(numout,*) ' : ',tname(jh),' ',ana_freq(jh)
END DO
! Initialize temporary arrays:
! ----------------------------
# if defined key_3Ddiaharm
ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 5, jpk ) )
ana_temp(:,:,:,:,:) = 0._wp
# else
ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 3 ) )
ana_temp(:,:,:,: ) = 0._wp
#endif
END SUBROUTINE dia_harm_init
SUBROUTINE dia_harm ( kt )
!!----------------------------------------------------------------------
!! *** ROUTINE dia_harm ***
!!
!! ** Purpose : Tidal harmonic analysis main routine
!!
!! ** Action : Sums ssh/u/v over time analysis [nit000_han,nitend_han]
!!
!!--------------------------------------------------------------------
INTEGER, INTENT( IN ) :: kt
!
INTEGER :: ji, jj, jh, jc, nhc
# if defined key_3Ddiaharm
INTEGER :: jk
# endif
REAL(wp) :: ztime, ztemp
!!--------------------------------------------------------------------
IF( nn_timing == 1 ) CALL timing_start('dia_harm')
IF( kt == nit000 ) CALL dia_harm_init
IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
ztime = (kt-nit000+1) * rdt
!IF(lwp) WRITE(numout,*) "ztime OLD", kt, ztime, sshn(25,25)
nhc = 0
DO jh = 1, nb_ana
DO jc = 1, 2
nhc = nhc+1
ztemp =( MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
& +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
! ssh, ub, vb are stored at the last level of 5d array
DO jj = 1,jpj
DO ji = 1,jpi
! Elevation and currents
# if defined key_3Ddiaharm
ana_temp(ji,jj,nhc,1,jpk) = ana_temp(ji,jj,nhc,1,jpk) + ztemp*sshn(ji,jj)*ssmask (ji,jj)
ana_temp(ji,jj,nhc,2,jpk) = ana_temp(ji,jj,nhc,2,jpk) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
ana_temp(ji,jj,nhc,3,jpk) = ana_temp(ji,jj,nhc,3,jpk) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
ana_temp(ji,jj,nhc,5,jpk) = ana_temp(ji,jj,nhc,5,jpk) &
& + ztemp*bfrva(ji,jj)*vn(ji,jj,mbkv(ji,jj))*ssvmask(ji,jj)
ana_temp(ji,jj,nhc,4,jpk) = ana_temp(ji,jj,nhc,4,jpk) &
& + ztemp*bfrua(ji,jj)*un(ji,jj,mbku(ji,jj))*ssumask(ji,jj)
# else
ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)
ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
# endif
END DO
END DO
!
# if defined key_3Ddiaharm
! 3d velocity and density:
DO jk=1,jpk-1
DO jj = 1,jpj
DO ji = 1,jpi
! density and velocity
ana_temp(ji,jj,nhc,1,jk) = ana_temp(ji,jj,nhc,1,jk) + ztemp*rhd(ji,jj,jk)
ana_temp(ji,jj,nhc,2,jk) = ana_temp(ji,jj,nhc,2,jk) + ztemp*(un(ji,jj,jk)-un_b(ji,jj)) &
& *umask(ji,jj,jk)
ana_temp(ji,jj,nhc,3,jk) = ana_temp(ji,jj,nhc,3,jk) + ztemp*(vn(ji,jj,jk)-vn_b(ji,jj)) &
& *vmask(ji,jj,jk)
ana_temp(ji,jj,nhc,4,jk) = ana_temp(ji,jj,nhc,4,jk) + ztemp*wn(ji,jj,jk)
ana_temp(ji,jj,nhc,5,jk) = ana_temp(ji,jj,nhc,5,jk) - 0.5*grav*ztemp*(rhd(ji,jj,jk)+rhd(ji,jj,jk+1) )/max(rn2(ji,jj,jk),1.e-8_wp)
! IF(jk<=mbathy(ji,jj) ) ana_temp(ji,jj,nhc,5,jk) = ana_temp(ji,jj,nhc,5,jk) - &
! & 0.5*grav*ztemp*(rhd(ji,jj,jk)+rhd(ji,jj,jk+1) )/max(rn2(ji,jj,jk),1.e-8_wp)
END DO
END DO
ENDDO
# endif
END DO
END DO
!
END IF
IF ( kt == nitend_han ) CALL dia_harm_end
IF( nn_timing == 1 ) CALL timing_stop('dia_harm')
END SUBROUTINE dia_harm
SUBROUTINE dia_harm_end
!!----------------------------------------------------------------------
!! *** ROUTINE diaharm_end ***
!!
!! ** Purpose : Compute the Real and Imaginary part of tidal constituents
!!
!! ** Action : Decompose the signal on the harmonic constituents
!!
!!--------------------------------------------------------------------
INTEGER :: ji, jj, jh, jc, jn, nhan, jl
# if defined key_3Ddiaharm
INTEGER :: jk
# endif
INTEGER :: ksp, kun, keq
REAL(wp) :: ztime, ztime_ini, ztime_end
REAL(wp) :: X1,X2
REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ana_amp
!!--------------------------------------------------------------------
CALL wrk_alloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis
ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis
nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
# if defined key_3Ddiaharm
ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana), &
& out_u (jpi,jpj,jpk,2*nb_ana), &
& out_v (jpi,jpj,jpk,2*nb_ana), &
& out_w (jpi,jpj,jpk,2*nb_ana), &
& out_dzi(jpi,jpj,jpk,2*nb_ana) )
# else
ALLOCATE( out_eta(jpi,jpj,2*nb_ana), &
& out_u (jpi,jpj,2*nb_ana), &
& out_v (jpi,jpj,2*nb_ana) )
# endif
IF(lwp) WRITE(numout,*) 'ANA F OLD', ft
IF(lwp) WRITE(numout,*) 'ANA U OLD', ut
IF(lwp) WRITE(numout,*) 'ANA V OLD', vt
ninco = 2*nb_ana
ksp = 0
keq = 0
DO jn = 1, nhan
ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1)
keq = keq + 1
kun = 0
DO jh = 1, nb_ana
DO jc = 1, 2
kun = kun + 1
ksp = ksp + 1
nisparse(ksp) = keq
njsparse(ksp) = kun
valuesparse(ksp) = ( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
& + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
END DO
END DO
END DO
nsparse = ksp
! Density and Elevation:
# if defined key_3Ddiaharm
DO jk=1,jpk
# endif
DO jj = 1, jpj
DO ji = 1, jpi
! Fill input array
kun = 0
DO jh = 1, nb_ana
DO jc = 1, 2
kun = kun + 1
# if defined key_3Ddiaharm
ztmp4(kun)=ana_temp(ji,jj,kun,1,jk)
# else
ztmp4(kun)=ana_temp(ji,jj,kun,1)
# endif
END DO
END DO
CALL SUR_DETERMINE(jj)
! Fill output array
DO jh = 1, nb_ana
ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
END DO
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
DO jh = 1, nb_ana
X1 = ana_amp(ji,jj,jh,1)
X2 =-ana_amp(ji,jj,jh,2)
# if defined key_3Ddiaharm
out_eta(ji,jj,jk,jh ) = X1 * tmask_i(ji,jj)
out_eta(ji,jj,jk,jh+nb_ana) = X2 * tmask_i(ji,jj)
# else
out_eta(ji,jj ,jh ) = X1 * tmask_i(ji,jj)
out_eta(ji,jj ,jh+nb_ana) = X2 * tmask_i(ji,jj)
# endif
END DO
END DO
END DO
! u-component of velocity
DO jj = 1, jpj
DO ji = 1, jpi
! Fill input array
kun=0
DO jh = 1,nb_ana
DO jc = 1,2
kun = kun + 1
# if defined key_3Ddiaharm
ztmp4(kun)=ana_temp(ji,jj,kun,2,jk)
# else
ztmp4(kun)=ana_temp(ji,jj,kun,2)
# endif
END DO
END DO
CALL SUR_DETERMINE(jj+1)
! Fill output array
DO jh = 1, nb_ana
ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
END DO
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
DO jh = 1, nb_ana
X1= ana_amp(ji,jj,jh,1)
X2=-ana_amp(ji,jj,jh,2)
# if defined key_3Ddiaharm
out_u(ji,jj,jk, jh) = X1 * ssumask(ji,jj)
out_u(ji,jj,jk,nb_ana+jh) = X2 * ssumask(ji,jj)
# else
out_u(ji,jj, jh) = X1 * ssumask(ji,jj)
out_u(ji,jj, nb_ana+jh) = X2 * ssumask(ji,jj)
# endif
ENDDO
ENDDO
ENDDO
! v- velocity
DO jj = 1, jpj
DO ji = 1, jpi
! Fill input array
kun=0
DO jh = 1,nb_ana
DO jc = 1,2
kun = kun + 1
# if defined key_3Ddiaharm
ztmp4(kun)=ana_temp(ji,jj,kun,3,jk)
# else
ztmp4(kun)=ana_temp(ji,jj,kun,3)
# endif
END DO
END DO
CALL SUR_DETERMINE(jj+1)
! Fill output array
DO jh = 1, nb_ana
ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
END DO
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
DO jh = 1, nb_ana
X1=ana_amp(ji,jj,jh,1)
X2=-ana_amp(ji,jj,jh,2)
# if defined key_3Ddiaharm
out_v(ji,jj,jk, jh)=X1 * ssvmask(ji,jj)
out_v(ji,jj,jk,nb_ana+jh)=X2 * ssvmask(ji,jj)
# else
out_v(ji,jj, jh)=X1 * ssvmask(ji,jj)
out_v(ji,jj, nb_ana+jh)=X2 * ssvmask(ji,jj)
# endif
END DO
END DO
END DO
# if defined key_3Ddiaharm
! w- velocity
DO jj = 1, jpj
DO ji = 1, jpi
! Fill input array
kun=0
DO jh = 1,nb_ana
DO jc = 1,2
kun = kun + 1
ztmp4(kun)=ana_temp(ji,jj,kun,4,jk)
END DO
END DO
CALL SUR_DETERMINE(jj+1)
! Fill output array
DO jh = 1, nb_ana
ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
END DO
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
DO jh = 1, nb_ana
X1=ana_amp(ji,jj,jh,1)
X2=-ana_amp(ji,jj,jh,2)
out_w(ji,jj,jk, jh)=X1 * tmask_i(ji,jj)
out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)
END DO
END DO
END DO
! dzi- isopycnal displacements
DO jj = 1, jpj
DO ji = 1, jpi
! Fill input array
kun=0
DO jh = 1,nb_ana
DO jc = 1,2
kun = kun + 1
ztmp4(kun)=ana_temp(ji,jj,kun,5,jk)
END DO
END DO
CALL SUR_DETERMINE(jj+1)
! Fill output array
DO jh = 1, nb_ana
ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
END DO
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
DO jh = 1, nb_ana
X1=ana_amp(ji,jj,jh,1)
X2=-ana_amp(ji,jj,jh,2)
out_dzi(ji,jj,jk, jh)=X1 * tmask_i(ji,jj)
out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)
END DO
END DO
END DO
ENDDO ! jk
# endif
CALL dia_wri_harm ! Write results in files
CALL wrk_dealloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
!
END SUBROUTINE dia_harm_end
SUBROUTINE dia_wri_harm
!!--------------------------------------------------------------------
!! *** ROUTINE dia_wri_harm ***
!!
!! ** Purpose : Write tidal harmonic analysis results in a netcdf file
!!--------------------------------------------------------------------
CHARACTER(LEN=lc) :: cltext
CHARACTER(LEN=lc) :: &
cdfile_name_T , & ! name of the file created (T-points)
cdfile_name_U , & ! name of the file created (U-points)
cdfile_name_V ! name of the file created (V-points)
INTEGER :: jh
# if defined key_3Ddiaharm
CHARACTER(LEN=lc) :: cdfile_name_W ! name of the file created (W-points)
INTEGER :: jk
REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im
REAL(WP), ALLOCATABLE, DIMENSION (:,:) :: z2real, z2im
# endif
!!----------------------------------------------------------------------
#if defined key_dimgout
cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'
cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'
cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'
# if defined key_3Ddiaharm
cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc'
# endif
#endif
IF(lwp) WRITE(numout,*) ' '
IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
#if defined key_dimgout
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Output files: ', TRIM(cdfile_name_T)
IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_U)
IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_V)
# if defined key_3Ddiaharm
IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_W)
# endif
#endif
IF(lwp) WRITE(numout,*) ' '
# if defined key_3Ddiaharm
ALLOCATE( z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj))
# endif
! A) density and elevation
!/////////////
!
#if defined key_dimgout
cltext='density amplitude and phase; elevation is level=jpk '
CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')
#else
# if defined key_3Ddiaharm
z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
# endif
DO jh = 1, nb_ana
# if defined key_3Ddiaharm
DO jk=1,jpkm1
z3real(:,:,jk)=out_eta(:,:,jk,jh)
z3im (:,:,jk)=out_eta(:,:,jk,jh+nb_ana)
ENDDO
z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana)
CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) )
CALL iom_put( TRIM(tname(jh))//'y_ro', z3im (:,:,:) )
CALL iom_put( TRIM(tname(jh))//'x' , z2real(:,: ) )
CALL iom_put( TRIM(tname(jh))//'y' , z2im (:,: ) )
# else
WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh))
CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
# endif
END DO
#endif
! B) u
!/////////
!
#if defined key_dimgout
cltext='3d u amplitude and phase; ubar is the last level'
CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')
#else
# if defined key_3Ddiaharm
z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
# endif
DO jh = 1, nb_ana
# if defined key_3Ddiaharm
DO jk=1,jpkm1
z3real(:,:,jk)=out_u(:,:,jk,jh)
z3im (:,:,jk)=out_u(:,:,jk,jh+nb_ana)
ENDDO
z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana)
CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) )
CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:) )
CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) )
CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:) )
z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana)
CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) )
CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:) )
# else
CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) )
CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) )
# endif
END DO
#endif
! C) v
!/////////
!
#if defined key_dimgout
cltext='3d v amplitude and phase; vbar is the last level'
CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')
#else
# if defined key_3Ddiaharm
z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
# endif
DO jh = 1, nb_ana
# if defined key_3Ddiaharm
DO jk=1,jpkm1
z3real(:,:,jk)=out_v(:,:,jk,jh)
z3im (:,:,jk)=out_v(:,:,jk,jh+nb_ana)
ENDDO
z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana)
CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) )
CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:) )
CALL iom_put( TRIM(tname(jh))//'x_v2d' , z2real(:,:) )
CALL iom_put( TRIM(tname(jh))//'y_v2d' , z2im (:,:) )
z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana)
CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) )
CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:) )
# else
CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) )
CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) )
# endif
END DO
#endif
! D) w
# if defined key_3Ddiaharm
# if defined key_dimgout
cltext='3d w amplitude and phase; vort_baro is the last level'
CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')
# else
DO jh = 1, nb_ana
DO jk=1,jpkm1
z3real(:,:,jk)=out_w(:,:,jk,jh)
z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana)
ENDDO
CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) )
CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) )
END DO
# endif
! E) dzi + tau_bot
# if defined key_dimgout
cltext='dzi=g*ro/N2 amplitude and phase'
CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')
# else
DO jh = 1, nb_ana
DO jk=1,jpkm1
z3real(:,:,jk)=out_dzi(:,:,jk,jh)
z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana)
ENDDO
CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) )
CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) )
END DO
# endif
# endif
!
# if defined key_3Ddiaharm
DEALLOCATE(z3real, z3im, z2real,z2im)
# endif
END SUBROUTINE dia_wri_harm
SUBROUTINE SUR_DETERMINE(init)
!!---------------------------------------------------------------------------------
!! *** ROUTINE SUR_DETERMINE ***
!!
!!
!!
!!---------------------------------------------------------------------------------
INTEGER, INTENT(in) :: init
!
INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
REAL(wp) :: zval1, zval2, zx1
REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
!---------------------------------------------------------------------------------
CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
CALL wrk_alloc( jpincomax , ipos2 , ipivot )
IF( init == 1 ) THEN
IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
!
ztmp3(:,:) = 0._wp
!
DO jk1_sd = 1, nsparse
DO jk2_sd = 1, nsparse
nisparse(jk2_sd) = nisparse(jk2_sd)
njsparse(jk2_sd) = njsparse(jk2_sd)
IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) &
& + valuesparse(jk1_sd)*valuesparse(jk2_sd)
ENDIF
END DO
END DO
!
DO jj_sd = 1 ,ninco
ipos1(jj_sd) = jj_sd
ipos2(jj_sd) = jj_sd
END DO
!
DO ji_sd = 1 , ninco
!
!find greatest non-zero pivot:
zval1 = ABS(ztmp3(ji_sd,ji_sd))
!
ipivot(ji_sd) = ji_sd
DO jj_sd = ji_sd, ninco
zval2 = ABS(ztmp3(ji_sd,jj_sd))
IF( zval2.GE.zval1 )THEN
ipivot(ji_sd) = jj_sd
zval1 = zval2
ENDIF
END DO
!
DO ji1_sd = 1, ninco
zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd)
zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd))
ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd)
ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
END DO
!
ipos2(ji_sd) = ipos1(ipivot(ji_sd))
ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
ipos1(ji_sd) = ipos2(ji_sd)
ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
zpivot(ji_sd) = ztmp3(ji_sd,ji_sd)
DO jj_sd = 1, ninco
ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
END DO
!
DO ji2_sd = ji_sd+1, ninco
zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
DO jj_sd=1,ninco
ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
END DO
END DO
!
END DO
!
ENDIF ! End init==1
DO ji_sd = 1, ninco
ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
DO ji2_sd = ji_sd+1, ninco
ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
END DO
END DO
!system solving:
ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
ji_sd = ninco
DO ji_sd = ninco-1, 1, -1
zx1 = 0._wp
DO jj_sd = ji_sd+1, ninco
zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
END DO
ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
END DO
DO jj_sd =1, ninco
ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
END DO
CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
CALL wrk_dealloc( jpincomax , ipos2 , ipivot )
!
END SUBROUTINE SUR_DETERMINE
#else
!!----------------------------------------------------------------------
!! Default case : Empty module
!!----------------------------------------------------------------------
LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .FALSE.
CONTAINS
SUBROUTINE dia_harm ( kt ) ! Empty routine
INTEGER, INTENT( IN ) :: kt
WRITE(*,*) 'dia_harm: you should not have seen this print'
END SUBROUTINE dia_harm
#endif
!!======================================================================
END MODULE diaharm
NAMELISTS_AND_FORTRAN_FILES/f_files/diaharm_fast.F90 0 → 100755
View file @ 868014cc
MODULE diaharm_fast
!!======================================================================
!! *** MODULE example ***
!! Ocean physics: On line harmonic analyser
!!
!!=====================================================================
#if defined key_diaharm_fast
!!----------------------------------------------------------------------
!! 'key_harm_ana' : Calculate harmonic analysis
!!----------------------------------------------------------------------
!! harm_ana :
!! harm_ana_init :
!! NB: 2017-12 : add 3D harmonic analysis of velocities
!! integration of Maria Luneva's development
!! 'key_3Ddiaharm'
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE iom
USE in_out_manager ! I/O units
USE phycst ! physical constants
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE bdy_oce ! ocean open boundary conditions
USE bdytides ! tidal bdy forcing
USE daymod ! calendar
USE tideini
USE restart
USE ioipsl, ONLY : ju2ymds ! for calendar
!
!
USE timing ! preformance summary
USE zdf_oce
IMPLICIT NONE
PRIVATE
!! * Routine accessibility
PUBLIC dia_harm_fast ! routine called in step.F90 module
LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast = .TRUE. ! to be run or not
LOGICAL, PUBLIC :: lk_diaharm_2D ! = .TRUE. ! to run 2d
LOGICAL, PUBLIC :: lk_diaharm_3D ! = .TRUE. ! to run 3d
!! * Module variables
INTEGER, PARAMETER :: nharm_max = jpmax_harmo ! max number of harmonics to be analysed
INTEGER, PARAMETER :: nhm_max = 2*nharm_max+1
INTEGER, PARAMETER :: nvab = 2 ! number of 3D variables
INTEGER :: nharm
INTEGER :: nhm
INTEGER :: & !!! ** toto namelist (namtoto) **
nflag = 1 ! default value of nflag
REAL(wp), DIMENSION(nharm_max) :: &
om_tide ! tidal frequencies ( rads/sec)
REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:) :: &
bzz,c,x ! work arrays
REAL(wp) :: cca,ssa,zm,bt,dd_cumul
!
REAL(wp), PUBLIC :: fjulday_startharm !: Julian Day since start of harmonic analysis
REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf ! nodel/phase corrections used by diaharmana
REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a
!
INTEGER :: nvar_2d, nvar_3d !: number of 2d and 3d variables to analyse
INTEGER, ALLOCATABLE,DIMENSION(:) :: m_posi_2d, m_posi_3d
! Name of variables used in the restart
CHARACTER( LEN = 10 ), DIMENSION(5), PARAMETER :: m_varName2d = (/'ssh','u2d','v2d','ubfr','vbfr'/)
CHARACTER( LEN = 10 ), DIMENSION(4), PARAMETER :: m_varName3d = (/'rho','u3d','v3d','w3d'/)
!
REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,: ) :: g_cosamp2D, g_sinamp2D, g_cumul_var2D
REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:) :: g_cosamp3D, g_sinamp3D, g_cumul_var3D
!
REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) :: g_out2D,h_out2D ! arrays for output
REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: g_out3D,h_out3D ! arrays for 3D output
!
! NAMELIST
LOGICAL, PUBLIC :: ln_diaharm_store !: =T Stores data for harmonic Analysis
LOGICAL, PUBLIC :: ln_diaharm_compute !: =T Compute harmonic Analysis
LOGICAL, PUBLIC :: ln_diaharm_read_restart !: =T Read restart from a previous run
LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d
INTEGER :: nb_ana ! Number of harmonics to analyse
CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_all ! INDEX within the full set of constituents (tide.h90)
INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_sub ! INDEX within the subset of constituents pass in input
!! * Substitutions
!!----------------------------------------------------------------------
!! OPA 9.0 , LOCEAN-IPSL (2005)
!! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
!! or TOP 1.0 , LOCEAN-IPSL (2005)
!! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dia_harm_fast( kt )
!!----------------------------------------------------------------------
!! *** ROUTINE harm_ana ***
!!
!! ** Purpose : Harmonic analyser
!!
!! ** Method :
!!
!! ** Action : - first action (share memory array/varible modified
!! in this routine
!! - second action .....
!! - .....
!!
!! References :
!! Give references if exist otherwise suppress these lines
!!
!! History :
!! 9.0 ! 03-08 (Autor Names) Original code
!! ! 02-08 (Author names) brief description of modifications
!!----------------------------------------------------------------------
!! * Modules used
!! * arguments
INTEGER, INTENT( in ) :: &
kt ! describe it!!!
!! * local declarations
INTEGER :: ji, jk, jj ! dummy loop arguments
INTEGER :: jh, i1, i2, jgrid
INTEGER :: j2d, j3d
REAL(WP) :: sec2start
!!--------------------------------------------------------------------
IF( nn_timing == 1 ) CALL timing_start( 'dia_harm_fast' )
IF( kt == nit000 ) CALL harm_ana_init ! Initialization (first time-step only)
IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
! this bit done every time step
nhm=2*nb_ana+1
c(1) = 1.0
sec2start = nint( (fjulday-fjulday_startharm)*86400._wp )
!IF(lwp) WRITE(numout,*) "ztime NEW", kt, sec2start, fjulday_startharm
DO jh=1,nb_ana
c(2*jh ) = anaf(jh)*cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
c(2*jh+1) = anaf(jh)*sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
ENDDO
!IF(lwp) WRITE(numout,*) "c init", c, "c end", sec2start, om_tide(1), anau(1), anav(1),"end nodal"
! CUMULATE
DO ji=1,jpi ! loop lon
DO jj=1,jpj ! loop lat
DO jh=1,nhm ! loop harmonic
DO j2d=1,nvar_2d
IF ( m_posi_2d(j2d) .eq. 1 ) dd_cumul = c(jh) * sshn(ji,jj) * ssmask (ji,jj) ! analysis elevation
IF ( m_posi_2d(j2d) .eq. 2 ) dd_cumul = c(jh) * un_b(ji,jj) * ssumask(ji,jj) ! analysis depth average velocities
IF ( m_posi_2d(j2d) .eq. 3 ) dd_cumul = c(jh) * vn_b(ji,jj) * ssvmask(ji,jj)
IF ( m_posi_2d(j2d) .eq. 4 ) dd_cumul = c(jh) * bfrua(ji,jj) * un(ji,jj,mbku(ji,jj)) * ssumask(ji,jj) ! analysis bottom friction
IF ( m_posi_2d(j2d) .eq. 5 ) dd_cumul = c(jh) * bfrva(ji,jj) * vn(ji,jj,mbkv(ji,jj)) * ssvmask(ji,jj)
g_cumul_var2D(jh,ji,jj,j2d) = g_cumul_var2D(jh,ji,jj,j2d) + dd_cumul
ENDDO
DO j3d=1,nvar_3d
DO jk=1,jpkm1
IF ( m_posi_3d(j3d) .eq. 1 ) dd_cumul = c(jh) * rhd(ji,jj,jk) * tmask(ji,jj,jk)
IF ( m_posi_3d(j3d) .eq. 2 ) dd_cumul = c(jh) * ( un(ji,jj,jk)-un_b(ji,jj) ) * umask(ji,jj,jk)
IF ( m_posi_3d(j3d) .eq. 3 ) dd_cumul = c(jh) * ( vn(ji,jj,jk)-vn_b(ji,jj) ) * vmask(ji,jj,jk)
IF ( m_posi_3d(j3d) .eq. 4 ) dd_cumul = c(jh) * wn(ji,jj,jk) * wmask(ji,jj,jk)
g_cumul_var3D(jh,ji,jj,jk,j3d) = g_cumul_var3D(jh,ji,jj,jk,j3d) + dd_cumul
ENDDO
ENDDO
ENDDO ! end loop harmonic
ENDDO ! end loop lat
ENDDO ! end loop lon
! Compute nodal factor cumulative cross-product
DO i1=1,nhm
DO i2=1,nhm
cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
ENDDO
ENDDO
! Output RESTART
IF( kt == nitrst ) THEN
CALL harm_rst_write(kt) ! Dump out data for a restarted run
ENDIF
! At End of run
IF ( kt == nitend ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
IF(lwp) WRITE(numout,*) '~~~~~~~~~'
IF( ln_diaharm_compute ) THEN
! INITIALISE TABLE TO 0
IF ( nvar_2d .gt. 0 ) THEN
g_cosamp2D = 0.0_wp
g_sinamp2D = 0.0_wp
ENDIF
IF ( nvar_3d .gt. 0 ) THEN
g_cosamp3D = 0.0_wp
g_sinamp3D = 0.0_wp
ENDIF
! FIRST OUTPUT 2D VARIABLES
DO jgrid=1,nvar_2d ! loop number of 2d variables (ssh, U2d, V2d, UVfric) to analyse harmonically
DO ji=1,jpi ! loop lon
DO jj=1,jpj ! loop lat
bt = 1.0_wp; bzz(:) = 0.0_wp
DO jh=1,nhm ! loop harmonic
bzz(jh) = g_cumul_var2D(jh,ji,jj,jgrid)
bt = bt*bzz(jh)
ENDDO
! Copy back original cumulated nodal factor
a(:,:) = cc(:,:)
! now do gaussian elimination of the system
! a * x = b
! the matrix x is (a0,a1,b1,a2,b2 ...)
! the matrix a and rhs b solved here for x
x=0.0_wp
IF(bt.ne.0.) THEN
CALL gelim( a, bzz, x, nhm )
! Backup output in variables
DO jh=1,nb_ana
g_cosamp2D(jh,ji,jj,jgrid) = x(jh*2 )
g_sinamp2D(jh,ji,jj,jgrid) = x(jh*2+1)
ENDDO
g_cosamp2D( 0,ji,jj,jgrid) = x(1)
g_sinamp2D( 0,ji,jj,jgrid) = 0.0_wp
ENDIF ! bt.ne.0.
ENDDO ! jj
ENDDO ! ji
ENDDO ! jgrid
! SECOND OUTPUT 3D VARIABLES
DO jgrid=1,nvar_3d ! loop number of 3d variables rho, U, V, W
DO jk=1,jpkm1 ! loop over vertical level
DO ji=1,jpi ! loop over lon
DO jj=1,jpj ! loop over lat
bt = 1.0_wp; bzz(:) = 0.0_wp
DO jh=1,nhm
bzz(jh) = g_cumul_var3D(jh,ji,jj,jk,jgrid)
bt = bt*bzz(jh)
ENDDO
! Copy back original cumulated nodal factor
a(:,:) = cc(:,:)
! now do gaussian elimination of the system
! a * x = b
! the matrix x is (a0,a1,b1,a2,b2 ...)
! the matrix a and rhs b solved here for x
x=0.0_wp
IF(bt.ne.0.) THEN
CALL gelim( a, bzz, x, nhm )
! Backup output in variables
DO jh=1,nb_ana
g_cosamp3D(jh,ji,jj,jk,jgrid) = x(jh*2 )
g_sinamp3D(jh,ji,jj,jk,jgrid) = x(jh*2+1)
ENDDO
g_cosamp3D ( 0,ji,jj,jk,jgrid) = x(1)
g_sinamp3D ( 0,ji,jj,jk,jgrid) = 0.0_wp
ENDIF ! bt.ne.0.
ENDDO ! jj
ENDDO ! ji
ENDDO ! jk
ENDDO ! jgrid
CALL harm_ana_out ! output analysis (last time step)
ELSE ! ln_harmana_compute = False
IF(lwp) WRITE(numout,*) " Skipping Computing harmonics at last step"
ENDIF ! ln_harmana_compute
ENDIF ! kt == nitend
ENDIF
IF( nn_timing == 1 ) CALL timing_stop( 'dia_harm_fast' )
END SUBROUTINE dia_harm_fast
SUBROUTINE harm_ana_init
!!----------------------------------------------------------------------
!! *** ROUTINE harm_ana_init ***
!!
!! ** Purpose : initialization of ....
!!
!! ** Method : blah blah blah ...
!!
!! ** input : Namlist namexa
!!
!! ** Action : ...
!!
!! history :
!! 9.0 ! 03-08 (Autor Names) Original code
!!----------------------------------------------------------------------
!! * local declarations
INTEGER :: ji, jk, jh ! dummy loop indices
INTEGER :: ios ! Local integer output status for namelist read
INTEGER :: k2d, k3d ! dummy number of analysis
NAMELIST/nam_diaharm_fast/ ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, tname
!!----------------------------------------------------------------------
lk_diaharm_2D = .TRUE. ! to run 2d
lk_diaharm_3D = .TRUE. ! to run 3d
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
IF(lwp) WRITE(numout,*) '~~~~~~~~~'
! GET NAMELIST DETAILS
REWIND( numnam_ref ) ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis
READ ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist', lwp )
REWIND( numnam_cfg ) ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis
READ ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist', lwp )
IF(lwm) WRITE ( numond, nam_diaharm_fast )
! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90
nb_ana = 0
DO jk=1,jpmax_harmo
DO ji=1,nb_harmo
IF(TRIM(tname(jk)) == Wave( ntide(ji) )%cname_tide ) THEN
nb_ana=nb_ana+1
ENDIF
END DO
END DO
!
IF(lwp) THEN
WRITE(numout,*) ' Namelist nam_diaharm_fast'
WRITE(numout,*) ' nb_ana = ', nb_ana
CALL flush(numout)
ENDIF
!
IF (nb_ana > nharm_max) THEN
IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop'
IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max
nstop = nstop + 1
ENDIF
ALLOCATE(ntide_all(nb_ana))
ALLOCATE(ntide_sub(nb_ana))
DO jk=1,nb_ana
DO ji=1,nb_harmo
IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN
ntide_sub(jk) = ji
ntide_all(jk) = ntide(ji)
EXIT
END IF
END DO
END DO
! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS
nvar_2d = 0; nvar_3d = 0
IF ( ln_ana_ssh ) nvar_2d = nvar_2d + 1 ! analysis elevation
IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2 ! analysis depth-averaged velocity
IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2 ! analysis bottom friction
IF ( ln_ana_rho ) nvar_3d = nvar_3d + 1 ! analysis density
IF ( ln_ana_uv3d ) nvar_3d = nvar_3d + 2 ! analysis 3D horizontal velocities
IF ( ln_ana_w3d ) nvar_3d = nvar_3d + 1 ! analysis 3D vertical velocity
! CHECK IF SOMETHING TO RUN
IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE. ! no 2d to run
IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE. ! no 3d to run
! IF ( nvar_2d .gt. 0 .and. nvar_3d .gt. 0 ) lk_diaharm_fast = .FALSE.
! IF ( .NOT. ln_diaharm_store ) lk_diaharm_fast = .FALSE.
IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
! DO ALLOCATIONS
IF ( lk_diaharm_2D ) THEN
ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj, nvar_2d) )
ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) )
ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) )
ALLOCATE( g_out2D (jpi,jpj) )
ALLOCATE( h_out2D (jpi,jpj) )
ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0
ENDIF
IF ( lk_diaharm_3D ) THEN
ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
ALLOCATE( g_out3D (jpi,jpj,jpk) )
ALLOCATE( h_out3D (jpi,jpj,jpk) )
ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0
ENDIF
ALLOCATE( cc(nb_ana*2+1,nb_ana*2+1) )
ALLOCATE( a (nb_ana*2+1,nb_ana*2+1) )
ALLOCATE( bzz(nb_ana*2+1) )
ALLOCATE( x (nb_ana*2+1) )
ALLOCATE( c (nb_ana*2+1) )
ALLOCATE( anau(nb_ana) )
ALLOCATE( anav(nb_ana) )
ALLOCATE( anaf(nb_ana) )
! END ALLOCATE
! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL
! MAKES THINGS EASIER AND FASTER LATER
! !!! UGLY !!!
jh = 1; k2d = 0;
IF ( ln_ana_ssh ) THEN
k2d = k2d + 1; m_posi_2d(k2d) = jh
IF(lwp) WRITE(numout,*) " - ssh harmonic analysis activated (ln_ana_ssh)"
ENDIF
jh = jh + 1
IF ( ln_ana_uvbar ) THEN
k2d = k2d + 1; m_posi_2d(k2d) = jh
jh = jh + 1
k2d = k2d + 1; m_posi_2d(k2d) = jh
IF(lwp) WRITE(numout,*) " - barotropic currents harmonic analysis activated (ln_ana_uvbar)"
ELSE
jh = jh + 1
ENDIF
jh = jh + 1
IF ( ln_ana_bfric ) THEN
k2d = k2d + 1; m_posi_2d(k2d) = jh
jh = jh + 1;
k2d = k2d + 1; m_posi_2d(k2d) = jh
IF(lwp) WRITE(numout,*) " - bottom friction harmonic analysis activated (ln_ana_vbfr)"
ELSE
jh = jh + 1
ENDIF
! and for 3D
jh = 1; k3d = 0;
IF ( ln_ana_rho ) THEN
k3d = k3d + 1; m_posi_3d(k3d) = jh
IF(lwp) WRITE(numout,*) " - 3D density harmonic analysis activated (ln_ana_rho)"
ENDIF
jh = jh + 1
IF ( ln_ana_uv3d ) THEN
k3d = k3d + 1; m_posi_3d(k3d) = jh
jh = jh + 1
k3d = k3d + 1; m_posi_3d(k3d) = jh
IF(lwp) WRITE(numout,*) " - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)"
ELSE
jh = jh + 1
ENDIF
jh = jh + 1
IF ( ln_ana_w3d ) THEN
k3d = k3d + 1; m_posi_3d(k3d) = jh
IF(lwp) WRITE(numout,*) " - 3D vertical currents harmonic analysis activated (ln_ana_w3d)"
ENDIF
! SELECT AND STORE FREQUENCIES
IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
DO jh=1,nb_ana
om_tide(jh) = omega_tide( ntide_sub(jh) )
IF(lwp) WRITE(numout,*) ' - ',tname(jh),' ',om_tide(jh)
ENDDO
! READ RESTART IF
IF ( ln_diaharm_read_restart ) THEN
IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run"
! Need to read in bssh bz, cc anau anav and anaf
call harm_rst_read ! This reads in from the previous day
! Currrently the data in in assci format
ELSE
IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh "
IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,: ) = 0.0_wp
IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp
cc = 0.0_wp
a (:,:) = 0.0_wp ! NB
bzz (:) = 0.0_wp
x (:) = 0.0_wp
c (:) = 0.0_wp
anau (:) = 0.0_wp
anav (:) = 0.0_wp
anaf (:) = 0.0_wp
DO jh = 1, nb_ana
anau(jh) = utide ( ntide_sub(jh) )
anav(jh) = v0tide( ntide_sub(jh) )
anaf(jh) = ftide ( ntide_sub(jh) )
END DO
fjulday_startharm=fjulday !Set this at very start and store
IF (lwp) THEN
WRITE(numout,*) '--------------------------'
WRITE(numout,*) ' - Output anaf for check'
WRITE(numout,*) 'ANA F', anaf
WRITE(numout,*) 'ANA U', anau
WRITE(numout,*) 'ANA V', anav
WRITE(numout,*) fjulday_startharm
WRITE(numout,*) '--------------------------'
ENDIF
ENDIF
ELSE
IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis"
ENDIF
END SUBROUTINE harm_ana_init
!
SUBROUTINE gelim (a,b,x,n)
!!----------------------------------------------------------------------
!! *** ROUTINE harm_ana ***
!!
!! ** Purpose : Guassian elimination
!!
!!
!! ** Action : - first action (share memory array/varible modified
!! in this routine
!! - second action .....
!! - .....
!!
!! References :
!! Give references if exist otherwise suppress these lines
!!
!! History :
implicit none
!
integer :: n
REAL(WP) :: b(nb_ana*2+1), a(nb_ana*2+1,nb_ana*2+1)
REAL(WP) :: x(nb_ana*2+1)
INTEGER :: row,col,prow,pivrow,rrow
REAL(WP) :: atemp
REAL(WP) :: pivot
REAL(WP) :: m
do row=1,n-1
pivrow=row
pivot=a(row,n-row+1)
do prow=row+1,n
if (abs(a(prow,n-row+1)).gt.abs(pivot) ) then
pivot=a(prow,n-row+1)
pivrow=prow
endif
enddo
! swap row and prow
if ( pivrow .ne. row ) then
atemp=b(pivrow)
b(pivrow)=b(row)
b(row)=atemp
do col=1,n
atemp=a(pivrow,col)
a(pivrow,col)=a(row,col)
a(row,col)=atemp
enddo
endif
do rrow=row+1,n
if (a(row,row).ne.0) then
m=-a(rrow,n-row+1)/a(row,n-row+1)
do col=1,n
a(rrow,col)=m*a(row,col)+a(rrow,col)
enddo
b(rrow)=m*b(row)+b(rrow)
endif
enddo
enddo
! back substitution now
x(1)=b(n)/a(n,1)
do row=n-1,1,-1
x(n-row+1)=b(row)
do col=1,(n-row)
x(n-row+1)=(x(n-row+1)-a(row,col)*x(col))
enddo
x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
enddo
return
END SUBROUTINE gelim
SUBROUTINE harm_ana_out
!!----------------------------------------------------------------------
!! *** ROUTINE harm_ana_init ***
!!
!! ** Purpose : initialization of ....
!!
!! ** Method : blah blah blah ...
!!
!! ** input : Namlist namexa
!!
!! ** Action : ...
!!
!! history :
!! 9.0 ! 03-08 (Autor Names) Original code
!!----------------------------------------------------------------------
USE dianam ! build name of file (routine)
!! * local declarations
INTEGER :: ji, jj, jk, jgrid, jh ! dummy loop indices
! INTEGER :: nh_T
! INTEGER :: nid_harm
! CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names
! CHARACTER (len=40) :: clhstnamu, clhstnamv ! temporary names
CHARACTER (len=40) :: suffix
! REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi ! temporary scalars
do jgrid=1,nvar_2d
do jh=1,nb_ana
h_out2D = 0.0
g_out2D = 0.0
do jj=1,nlcj
do ji=1,nlci
cca=g_cosamp2D(jh,ji,jj,jgrid)
ssa=g_sinamp2D(jh,ji,jj,jgrid)
h_out2D(ji,jj)=sqrt(cca**2+ssa**2)
IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
g_out2D(ji,jj)= 0.0_wp
ELSE
g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca)
ENDIF
IF (h_out2D(ji,jj).ne.0) THEN
h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh)
ENDIF
IF (g_out2D(ji,jj).ne.0) THEN !Correct and take modulus
g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
if (g_out2D(ji,jj).gt.360.0) then
g_out2D(ji,jj)=g_out2D(ji,jj)-360.0
else if (g_out2D(ji,jj).lt.0.0) then
g_out2D(ji,jj)=g_out2D(ji,jj)+360.0
endif
ENDIF
enddo
enddo
!
! NETCDF OUTPUT
suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix), h_out2D(:,:) )
CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix), g_out2D(:,:) )
enddo
enddo
!
! DO THE SAME FOR 3D VARIABLES
!
do jgrid=1,nvar_3d
do jh=1,nb_ana
h_out3D = 0.0
g_out3D = 0.0
DO jk=1,jpkm1
do jj=1,nlcj
do ji=1,nlci
cca=g_cosamp3D(jh,ji,jj,jk,jgrid)
ssa=g_sinamp3D(jh,ji,jj,jk,jgrid)
h_out3D(ji,jj,jk)=sqrt(cca**2+ssa**2)
IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
g_out3D(ji,jj,jk) = 0.0_wp
ELSE
g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca)
ENDIF
IF (h_out3D(ji,jj,jk).ne.0) THEN
h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh)
ENDIF
IF (g_out3D(ji,jj,jk).ne.0) THEN !Correct and take modulus
g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
if (g_out3D(ji,jj,jk).gt.360.0) then
g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0
else if (g_out3D(ji,jj,jk).lt.0.0) then
g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0
endif
ENDIF
enddo ! ji
enddo ! jj
ENDDO ! jk
!
! NETCDF OUTPUT
suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix), h_out3D(:,:,:) )
CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix), g_out3D(:,:,:) )
enddo ! jh
enddo ! jgrid
!
END SUBROUTINE harm_ana_out
!
SUBROUTINE harm_rst_write(kt)
!!----------------------------------------------------------------------
!! *** ROUTINE harm_ana_init ***
!!
!! ** Purpose : To write out cummulated Tidal Harmomnic data to file for
!! restarting
!!
!! ** Method : restart files will be dated by default
!!
!! ** input :
!!
!! ** Action : ...
!!
!! history :
!! 0.0 ! 01-16 (Enda O'Dea) Original code
!! ASSUMES dated file for rose , can change later to be more generic
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step
!!
INTEGER :: jh, j2d, j3d
CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character
CHARACTER(LEN=50) :: clname ! ocean output restart file name
CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file
CHARACTER(LEN=250) :: clfinal ! full name
!restart file
DO j2d=1,nvar_2d
CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
DO jh=1,nb_ana
CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2 , :, :, j2d ) )
CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
ENDDO
ENDDO
DO j3d=1,nvar_3d
CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
DO jh=1,nb_ana
CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2 , :, :, :, j3d ) )
CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
ENDDO
ENDDO
IF(lwp) THEN
IF( kt > 999999999 ) THEN ; WRITE(clkt, * ) kt
ELSE ; WRITE(clkt, '(i8.8)') kt
ENDIF
clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
clpath = TRIM(cn_ocerst_outdir)
IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
WRITE(66) cc
WRITE(66) anau
WRITE(66) anav
WRITE(66) anaf
WRITE(66) fjulday_startharm
CLOSE(66)
WRITE(numout,*) '----------------------------'
WRITE(numout,*) ' harm_rst_write: DONE '
WRITE(numout,*) cc
WRITE(numout,*) anaf
WRITE(numout,*) fjulday_startharm
WRITE(numout,*) '----------------------------'
ENDIF
END SUBROUTINE harm_rst_write
SUBROUTINE harm_rst_read
!!----------------------------------------------------------------------
!! *** ROUTINE harm_ana_init ***
!!
!! ** Purpose : To read in cummulated Tidal Harmomnic data to file for
!! restarting
!!
!! ** Method :
!!
!! ** input :
!!
!! ** Action : ...
!!
!! history :
!! 0.0 ! 01-16 (Enda O'Dea) Original code
!! ASSUMES dated file for rose , can change later to be more generic
!!----------------------------------------------------------------------
CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character
CHARACTER(LEN=50) :: clname ! ocean output restart file name
CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file
CHARACTER(LEN=250) :: clfinal ! full name
INTEGER :: jh, j2d, j3d
IF( nit000 > 999999999 ) THEN ; WRITE(clkt, * ) nit000-1
ELSE ; WRITE(clkt, '(i8.8)') nit000-1
ENDIF
clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
clpath = TRIM(cn_ocerst_outdir)
IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
DO j2d=1,nvar_2d
CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
IF(lwp) WRITE(numout,*) "2D", j2d, m_posi_2d(j2d), m_varName2d( m_posi_2d(j2d) )
DO jh=1,nb_ana
CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2 , :, :, j2d ) )
CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
ENDDO
ENDDO
DO j3d=1,nvar_3d
CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
IF(lwp) WRITE(numout,*) "3D", j3d, m_posi_3d(j3d), m_varName3d( m_posi_3d(j3d) )
DO jh=1,nb_ana
CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2 , :, :, :, j3d ) )
CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
ENDDO
ENDDO
WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
READ(66) cc
READ(66) anau
READ(66) anav
READ(66) anaf
READ(66) fjulday_startharm
CLOSE(66)
IF(lwp) THEN
WRITE(numout,*) '----------------------------'
WRITE(numout,*) ' Checking anaf is correct'
WRITE(numout,*) cc
WRITE(numout,*) anaf
WRITE(numout,*) fjulday_startharm
WRITE(numout,*) '----------------------------'
ENDIF
END SUBROUTINE harm_rst_read
!!======================================================================
#else
!!---------------------------------------------------------------------------------
!! Dummy module NO harmonic Analysis
!!---------------------------------------------------------------------------------
LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast = .FALSE. ! to be run or not
CONTAINS
SUBROUTINE harm_rst_write(kt) ! Dummy routine
END SUBROUTINE harm_rst_write
SUBROUTINE harm_rst_read ! Dummy routine
END SUBROUTINE harm_rst_read
SUBROUTINE harm_ana_out ! Dummy routine
END SUBROUTINE harm_ana_out
SUBROUTINE harm_ana_init
END SUBROUTINE harm_ana_init
SUBROUTINE harm_ana( kt )
!--- NB : end call not properly written
END SUBROUTINE harm_ana
! END SUBROUTINE harm_ana_init
!--- END NB
SUBROUTINE gelim (a,b,x,n)
!--- NB : end call not properly written
END SUBROUTINE gelim
! END SUBROUTINE gelim (a,b,x,n)
!--- END NB
#endif
END MODULE diaharm_fast
NAMELISTS_AND_FORTRAN_FILES/f_files/dommsk.F90 0 → 100755
View file @ 868014cc
MODULE dommsk
!!======================================================================
!! *** MODULE dommsk ***
!! Ocean initialization : domain land/sea mask
!!======================================================================
!! History : OPA ! 1987-07 (G. Madec) Original code
!! 6.0 ! 1993-03 (M. Guyon) symetrical conditions (M. Guyon)
!! 7.0 ! 1996-01 (G. Madec) suppression of common work arrays
!! - ! 1996-05 (G. Madec) mask computed from tmask
!! 8.0 ! 1997-02 (G. Madec) mesh information put in domhgr.F
!! 8.1 ! 1997-07 (G. Madec) modification of kbat and fmask
!! - ! 1998-05 (G. Roullet) free surface
!! 8.2 ! 2000-03 (G. Madec) no slip accurate
!! - ! 2001-09 (J.-M. Molines) Open boundaries
!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
!! - ! 2005-11 (V. Garnier) Surface pressure gradient organization
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!! 3.6 ! 2015-05 (P. Mathiot) ISF: add wmask,wumask and wvmask
!! 4.0 ! 2016-06 (G. Madec, S. Flavoni) domain configuration / user defined interface
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! dom_msk : compute land/ocean mask
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE usrdef_fmask ! user defined fmask
USE bdy_oce
USE in_out_manager ! I/O manager
USE iom
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp ! Massively Parallel Processing library
USE wrk_nemo ! Memory allocation
USE timing ! Timing
IMPLICIT NONE
PRIVATE
PUBLIC dom_msk ! routine called by inidom.F90
! !!* Namelist namlbc : lateral boundary condition *
REAL(wp) :: rn_shlat ! type of lateral boundary condition on velocity
LOGICAL, PUBLIC :: ln_vorlat ! consistency of vorticity boundary condition
! with analytical eqs.
!! * Substitutions
# include "vectopt_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OPA 3.2 , LODYC-IPSL (2009)
!! $Id: dommsk.F90 7753 2017-03-03 11:46:59Z mocavero $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dom_msk( k_top, k_bot )
!!---------------------------------------------------------------------
!! *** ROUTINE dom_msk ***
!!
!! ** Purpose : Compute land/ocean mask arrays at tracer points, hori-
!! zontal velocity points (u & v), vorticity points (f) points.
!!
!! ** Method : The ocean/land mask at t-point is deduced from ko_top
!! and ko_bot, the indices of the fist and last ocean t-levels which
!! are either defined in usrdef_zgr or read in zgr_read.
!! The velocity masks (umask, vmask, wmask, wumask, wvmask)
!! are deduced from a product of the two neighboring tmask.
!! The vorticity mask (fmask) is deduced from tmask taking
!! into account the choice of lateral boundary condition (rn_shlat) :
!! rn_shlat = 0, free slip (no shear along the coast)
!! rn_shlat = 2, no slip (specified zero velocity at the coast)
!! 0 < rn_shlat < 2, partial slip | non-linear velocity profile
!! 2 < rn_shlat, strong slip | in the lateral boundary layer
!!
!! tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
!! rows/lines due to cyclic or North Fold boundaries as well
!! as MPP halos.
!! tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines
!! due to cyclic or North Fold boundaries as well as MPP halos.
!!
!! ** Action : tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask
!! at t-, u-, v- w, wu-, and wv-points (=0. or 1.)
!! fmask : land/ocean mask at f-point (=0., or =1., or
!! =rn_shlat along lateral boundaries)
!! tmask_i : interior ocean mask
!! tmask_h : halo mask
!! ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask
!!----------------------------------------------------------------------
INTEGER, DIMENSION(:,:), INTENT(in) :: k_top, k_bot ! first and last ocean level
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: iif, iil ! local integers
INTEGER :: ijf, ijl ! - -
INTEGER :: iktop, ikbot ! - -
INTEGER :: ios, inum
REAL(wp), POINTER, DIMENSION(:,:) :: zwf ! 2D workspace
!!
NAMELIST/namlbc/ rn_shlat, ln_vorlat
NAMELIST/nambdy/ ln_bdy ,nb_bdy, ln_coords_file, cn_coords_file, &
& ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
& cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
& ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
& cn_ice_lim, nn_ice_lim_dta, &
& rn_ice_tem, rn_ice_sal, rn_ice_age, &
& ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
!!---------------------------------------------------------------------
!
IF( nn_timing == 1 ) CALL timing_start('dom_msk')
!
REWIND( numnam_ref ) ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
READ ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
REWIND( numnam_cfg ) ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
READ ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
IF(lwm) WRITE ( numond, namlbc )
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) 'dommsk : ocean mask '
WRITE(numout,*) '~~~~~~'
WRITE(numout,*) ' Namelist namlbc'
WRITE(numout,*) ' lateral momentum boundary cond. rn_shlat = ',rn_shlat
WRITE(numout,*) ' consistency with analytical form ln_vorlat = ',ln_vorlat
ENDIF
IF ( rn_shlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral free-slip '
ELSEIF ( rn_shlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral no-slip '
ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral partial-slip '
ELSEIF ( 2. < rn_shlat ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral strong-slip '
ELSE
WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
CALL ctl_stop( ctmp1 )
ENDIF
! Ocean/land mask at t-point (computed from ko_top and ko_bot)
! ----------------------------
!
tmask(:,:,:) = 0._wp
DO jj = 1, jpj
DO ji = 1, jpi
iktop = k_top(ji,jj)
ikbot = k_bot(ji,jj)
IF( iktop /= 0 ) THEN ! water in the column
tmask(ji,jj,iktop:ikbot ) = 1._wp
ENDIF
END DO
END DO
!SF add here lbc_lnk: bug not still understood : cause now domain configuration is read !
!!gm I don't understand why...
CALL lbc_lnk( tmask , 'T', 1._wp ) ! Lateral boundary conditions
! Mask corrections for bdy (read in mppini2)
REWIND( numnam_ref ) ! Namelist nambdy in reference namelist :Unstructured open boundaries
READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist :Unstructured open boundaries
READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
! ------------------------
IF ( ln_bdy .AND. ln_mask_file ) THEN
CALL iom_open( cn_mask_file, inum )
CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
CALL iom_close( inum )
DO jk = 1, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
END DO
END DO
END DO
ENDIF
! Ocean/land mask at u-, v-, and f-points (computed from tmask)
! ----------------------------------------
! NB: at this point, fmask is designed for free slip lateral boundary condition
DO jk = 1, jpk
DO jj = 1, jpjm1
DO ji = 1, fs_jpim1 ! vector loop
umask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk)
vmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji ,jj+1,jk)
END DO
DO ji = 1, jpim1 ! NO vector opt.
fmask(ji,jj,jk) = tmask(ji,jj ,jk) * tmask(ji+1,jj ,jk) &
& * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
END DO
END DO
END DO
CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions
CALL lbc_lnk( vmask , 'V', 1._wp )
CALL lbc_lnk( fmask , 'F', 1._wp )
! Ocean/land mask at wu-, wv- and w points (computed from tmask)
!-----------------------------------------
wmask (:,:,1) = tmask(:,:,1) ! surface
wumask(:,:,1) = umask(:,:,1)
wvmask(:,:,1) = vmask(:,:,1)
DO jk = 2, jpk ! interior values
wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)
wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
END DO
! Ocean/land column mask at t-, u-, and v-points (i.e. at least 1 wet cell in the vertical)
! ----------------------------------------------
ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 )
ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 )
ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
! Interior domain mask (used for global sum)
! --------------------
!
iif = jpreci ; iil = nlci - jpreci + 1
ijf = jprecj ; ijl = nlcj - jprecj + 1
!
! ! halo mask : 0 on the halo and 1 elsewhere
tmask_h(:,:) = 1._wp
tmask_h( 1 :iif, : ) = 0._wp ! first columns
tmask_h(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns)
tmask_h( : , 1 :ijf) = 0._wp ! first rows
tmask_h( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows)
!
! ! north fold mask
tpol(1:jpiglo) = 1._wp
fpol(1:jpiglo) = 1._wp
IF( jperio == 3 .OR. jperio == 4 ) THEN ! T-point pivot
tpol(jpiglo/2+1:jpiglo) = 0._wp
fpol( 1 :jpiglo) = 0._wp
IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row for tmask_h
DO ji = iif+1, iil-1
tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
END DO
ENDIF
ENDIF
!
IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot
tpol( 1 :jpiglo) = 0._wp
fpol(jpiglo/2+1:jpiglo) = 0._wp
ENDIF
!
! ! interior mask : 2D ocean mask x halo mask
tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:)
! Lateral boundary conditions on velocity (modify fmask)
! ---------------------------------------
IF( rn_shlat /= 0 ) THEN ! Not free-slip lateral boundary condition
!
CALL wrk_alloc( jpi,jpj, zwf )
!
DO jk = 1, jpk
zwf(:,:) = fmask(:,:,jk)
DO jj = 2, jpjm1
DO ji = fs_2, fs_jpim1 ! vector opt.
IF( fmask(ji,jj,jk) == 0._wp ) THEN
fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), &
& zwf(ji-1,jj), zwf(ji,jj-1) ) )
ENDIF
END DO
END DO
DO jj = 2, jpjm1
IF( fmask(1,jj,jk) == 0._wp ) THEN
fmask(1 ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
ENDIF
IF( fmask(jpi,jj,jk) == 0._wp ) THEN
fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
ENDIF
END DO
DO ji = 2, jpim1
IF( fmask(ji,1,jk) == 0._wp ) THEN
fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
ENDIF
IF( fmask(ji,jpj,jk) == 0._wp ) THEN
fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
ENDIF
END DO
END DO
!
CALL wrk_dealloc( jpi,jpj, zwf )
!
CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
!
! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) depending on ln_vorlat
!
ENDIF
! User defined alteration of fmask (use to reduce ocean transport in specified straits)
! --------------------------------
!
CALL usr_def_fmask( cn_cfg, nn_cfg, fmask )
!
!
IF( nn_timing == 1 ) CALL timing_stop('dom_msk')
!
END SUBROUTINE dom_msk
!!======================================================================
END MODULE dommsk
NAMELISTS_AND_FORTRAN_FILES/f_files/dtatsd.F90 0 → 100755
View file @ 868014cc
MODULE dtatsd
!!======================================================================
!! *** MODULE dtatsd ***
!! Ocean data : read ocean Temperature & Salinity Data from gridded data
!!======================================================================
!! History : OPA ! 1991-03 () Original code
!! - ! 1992-07 (M. Imbard)
!! 8.0 ! 1999-10 (M.A. Foujols, M. Imbard) NetCDF FORMAT
!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
!! 3.3 ! 2010-10 (C. Bricaud, S. Masson) use of fldread
!! 3.4 ! 2010-11 (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! dta_tsd : read and time interpolated ocean Temperature & Salinity Data
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE fldread ! read input fields
USE in_out_manager ! I/O manager
USE phycst ! physical constants
USE lib_mpp ! MPP library
USE wrk_nemo ! Memory allocation
USE timing ! Timing
USE iom
IMPLICIT NONE
PRIVATE
PUBLIC dta_tsd_init ! called by opa.F90
PUBLIC dta_tsd ! called by istate.F90 and tradmp.90
LOGICAL , PUBLIC :: ln_tsd_init !: T & S data flag
LOGICAL , PUBLIC :: ln_tsd_interp !: vertical interpolation flag
LOGICAL , PUBLIC :: ln_tsd_tradmp !: internal damping toward input data flag
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsd ! structure of input SST (file informations, fields read)
INTEGER :: jpk_init , inum_dta
INTEGER :: id ,linum ! local integers
INTEGER :: zdim(4)
!!----------------------------------------------------------------------
!! NEMO/OPA 3.3 , NEMO Consortium (2010)
!! $Id: dtatsd.F90 7753 2017-03-03 11:46:59Z mocavero $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dta_tsd_init( ld_tradmp )
!!----------------------------------------------------------------------
!! *** ROUTINE dta_tsd_init ***
!!
!! ** Purpose : initialisation of T & S input data
!!
!! ** Method : - Read namtsd namelist
!! - allocates T & S data structure
!!----------------------------------------------------------------------
LOGICAL, INTENT(in), OPTIONAL :: ld_tradmp ! force the initialization when tradp is used
!
INTEGER :: ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5 ! local integers
!!
CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
TYPE(FLD_N), DIMENSION(jpts+2):: slf_i ! array of namelist informations on the fields to read
TYPE(FLD_N) :: sn_tem, sn_sal, sn_dep, sn_msk
!!
NAMELIST/namtsd/ ln_tsd_init, ln_tsd_interp, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal, sn_dep, sn_msk
!!----------------------------------------------------------------------
!
IF( nn_timing == 1 ) CALL timing_start('dta_tsd_init')
!
! Initialisation
ierr0 = 0 ; ierr1 = 0 ; ierr2 = 0 ; ierr3 = 0 ; ierr4 = 0 ; ierr5 = 0
!
REWIND( numnam_ref ) ! Namelist namtsd in reference namelist :
READ ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in reference namelist', lwp )
REWIND( numnam_cfg ) ! Namelist namtsd in configuration namelist : Parameters of the run
READ ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp )
IF(lwm) WRITE ( numond, namtsd )
IF( PRESENT( ld_tradmp ) ) ln_tsd_tradmp = .TRUE. ! forces the initialization when tradmp is used
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) 'dta_tsd_init : Temperature & Salinity data '
WRITE(numout,*) '~~~~~~~~~~~~ '
WRITE(numout,*) ' Namelist namtsd'
WRITE(numout,*) ' Initialisation of ocean T & S with T &S input data ln_tsd_init = ', ln_tsd_init
WRITE(numout,*) ' iInterpolation of initial conditions in the vertical ln_tsd_interp = ', ln_tsd_interp
WRITE(numout,*) ' damping of ocean T & S toward T &S input data ln_tsd_tradmp = ', ln_tsd_tradmp
WRITE(numout,*)
IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_tradmp ) THEN
WRITE(numout,*)
WRITE(numout,*) ' T & S data not used'
ENDIF
ENDIF
!
IF( ln_rstart .AND. ln_tsd_init ) THEN
CALL ctl_warn( 'dta_tsd_init: ocean restart and T & S data intialisation, ', &
& 'we keep the restart T & S values and set ln_tsd_init to FALSE' )
ln_tsd_init = .FALSE.
ENDIF
IF( ln_tsd_interp .AND. ln_tsd_tradmp ) THEN
CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' ) ; RETURN
ENDIF
IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN
CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' ) ; RETURN
ENDIF
!
! ! allocate the arrays (if necessary)
IF( ln_tsd_init .OR. ln_tsd_tradmp ) THEN
!
IF( ln_tsd_interp ) THEN
ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information
ELSE
ALLOCATE( sf_tsd(jpts ), STAT=ierr0 )
ENDIF
IF( ierr0 > 0 ) THEN
CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' ) ; RETURN
ENDIF
!
IF( ln_tsd_interp ) THEN
CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta )
id = iom_varid( inum_dta, sn_dep%clvar, zdim )
jpk_init = zdim(3)
IF(lwp) WRITE(numout,*) 'Dimension of veritcal coordinate in ICs: ', jpk_init
CALL iom_close( inum_dta ) ! Close the input file
!
ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init ) , STAT=ierr0 )
IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 )
ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init ) , STAT=ierr2 )
IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 )
ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init ) , STAT=ierr4 )
ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init ) , STAT=ierr5 )
ELSE
ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 )
IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
ENDIF ! ln_tsd_interp
!
IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN
CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' ) ; RETURN
ENDIF
! ! fill sf_tsd with sn_tem & sn_sal and control print
slf_i(jp_tem) = sn_tem ; slf_i(jp_sal) = sn_sal
IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep ; slf_i(jp_msk) = sn_msk
CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print )
!
ENDIF
!
IF( nn_timing == 1 ) CALL timing_stop('dta_tsd_init')
!
END SUBROUTINE dta_tsd_init
SUBROUTINE dta_tsd( kt, ptsd )
!!----------------------------------------------------------------------
!! *** ROUTINE dta_tsd ***
!!
!! ** Purpose : provides T and S data at kt
!!
!! ** Method : - call fldread routine
!! - ORCA_R2: add some hand made alteration to read data
!! - 'key_orca_lev10' interpolates on 10 times more levels
!! - s- or mixed z-s coordinate: vertical interpolation on model mesh
!! - ln_tsd_tradmp=F: deallocates the T-S data structure
!! as T-S data are no are used
!!
!! ** Action : ptsd T-S data on medl mesh and interpolated at time-step kt
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: ptsd ! T & S data
!
INTEGER :: ji, jj, jk, jl, jk_init ! dummy loop indicies
INTEGER :: ik, il0, il1, ii0, ii1, ij0, ij1 ! local integers
REAL(wp):: zl, zi
!!----------------------------------------------------------------------
!
IF( nn_timing == 1 ) CALL timing_start('dta_tsd')
!
CALL fld_read( kt, 1, sf_tsd ) !== read T & S data at kt time step ==!
!
!
!!gm This should be removed from the code ===>>>> T & S files has to be changed
!
! !== ORCA_R2 configuration and T & S damping ==!
IF( cn_cfg == "orca" .AND. nn_cfg == 2 .AND. ln_tsd_tradmp ) THEN ! some hand made alterations
!
ij0 = 101 ; ij1 = 109 ! Reduced T & S in the Alboran Sea
ii0 = 141 ; ii1 = 155
DO jj = mj0(ij0), mj1(ij1)
DO ji = mi0(ii0), mi1(ii1)
sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
!
sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
sf_tsd(jp_sal)%fnow(ji,jj,18:25) = sf_tsd(jp_sal)%fnow(ji,jj,18:25) - 0.35_wp
END DO
END DO
ij0 = 87 ; ij1 = 96 ! Reduced temperature in Red Sea
ii0 = 148 ; ii1 = 160
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 4:10 ) = 7.0_wp
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
ENDIF
!!gm end
!
IF( kt == nit000 .AND. lwp )THEN
WRITE(numout,*)
WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh'
ENDIF
!
IF( ln_tsd_interp ) THEN ! probably should use pointers in the following to make more readable
!
DO jk = 1, jpk ! determines the intepolated T-S profiles at each (i,j) points
DO jj= 1, jpj
DO ji= 1, jpi
zl = gdept_0(ji,jj,jk)
IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN ! above the first level of data
ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1)
ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1)
ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN ! below the last level of data
ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init)
ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init)
ELSE ! inbetween : vertical interpolation between jk_init & jk_init+1
DO jk_init = 1, jpk_init-1 ! when gdept(jk_init) < zl < gdept(jk_init+1)
IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN ! if there is no data fill down
sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init)
sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init)
ENDIF
IF( (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN
zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / &
& (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init))
ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + &
& (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi
ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + &
& (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi
ENDIF
END DO
ENDIF
ENDDO
ENDDO
END DO
!
ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:)
ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:)
ELSE !== z- or zps- coordinate ==!
!
ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:) * tmask(:,:,:) ! Mask
ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) * tmask(:,:,:)
!
IF( ln_zps ) THEN ! zps-coordinate (partial steps) interpolation at the last ocean level
DO jj = 1, jpj
DO ji = 1, jpi
ik = mbkt(ji,jj)
IF( ik > 1 ) THEN
zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
ENDIF
ik = mikt(ji,jj)
IF( ik > 1 ) THEN
zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
END IF
END DO
END DO
ENDIF
!
ENDIF
!
IF( .NOT.ln_tsd_tradmp ) THEN !== deallocate T & S structure ==!
! (data used only for initialisation)
IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run'
DEALLOCATE( sf_tsd(jp_tem)%fnow ) ! T arrays in the structure
IF( sf_tsd(jp_tem)%ln_tint ) DEALLOCATE( sf_tsd(jp_tem)%fdta )
DEALLOCATE( sf_tsd(jp_sal)%fnow ) ! S arrays in the structure
IF( sf_tsd(jp_sal)%ln_tint ) DEALLOCATE( sf_tsd(jp_sal)%fdta )
IF( ln_tsd_interp ) DEALLOCATE( sf_tsd(jp_dep)%fnow ) ! T arrays in the structure
IF( ln_tsd_interp ) DEALLOCATE( sf_tsd(jp_msk)%fnow ) ! T arrays in the structure
DEALLOCATE( sf_tsd ) ! the structure itself
ENDIF
!
IF( nn_timing == 1 ) CALL timing_stop('dta_tsd')
!
END SUBROUTINE dta_tsd
!!======================================================================
END MODULE dtatsd
NAMELISTS_AND_FORTRAN_FILES/f_files/par_oce.F90 0 → 100755
View file @ 868014cc
MODULE par_oce
!!======================================================================
!! *** par_oce ***
!! Ocean : set the ocean parameters
!!======================================================================
!! History : OPA ! 1991 (Imbard, Levy, Madec) Original code
!! NEMO 1.0 ! 2004-01 (G. Madec, J.-M. Molines) Free form and module
!! 3.3 ! 2010-09 (C. Ethe) TRA-TRC merge: add jpts, jp_tem & jp_sal
!!----------------------------------------------------------------------
USE par_kind ! kind parameters
IMPLICIT NONE
PUBLIC
!!----------------------------------------------------------------------
!! namcfg namelist parameters
!!----------------------------------------------------------------------
LOGICAL :: ln_read_cfg !: (=T) read the domain configuration file or (=F) not
CHARACTER(lc) :: cn_domcfg !: filename the configuration file to be read
LOGICAL :: ln_write_cfg !: (=T) create the domain configuration file
CHARACTER(lc) :: cn_domcfg_out !: filename the configuration file to be read
!
LOGICAL :: ln_use_jattr !: input file read offset
! ! Use file global attribute: open_ocean_jstart to determine start j-row
! ! when reading input from those netcdf files that have the
! ! attribute defined. This is designed to enable input files associated
! ! with the extended grids used in the under ice shelf configurations to
! ! be used without redundant rows when the ice shelves are not in use.
!
!!---------------------------------------------------------------------
!! Domain Matrix size
!!---------------------------------------------------------------------
! configuration name & resolution (required only in ORCA family case)
CHARACTER(lc) :: cn_cfg !: name of the configuration
INTEGER :: nn_cfg !: resolution of the configuration
! global domain size !!! * total computational domain *
INTEGER :: jpiglo !: 1st dimension of global domain --> i-direction
INTEGER :: jpjglo !: 2nd - - --> j-direction
INTEGER :: jpkglo !: 3nd - - --> k levels
#if defined key_agrif
!!gm BUG ? I'm surprised by the calculation below of nbcellsx and nbcellsy before jpiglo,jpjglo
!!gm has been assigned to a value....
!!gm
! global domain size for AGRIF !!! * total AGRIF computational domain *
INTEGER, PUBLIC, PARAMETER :: nbghostcells = 1 !: number of ghost cells
INTEGER, PUBLIC :: nbcellsx = jpiglo - 2 - 2*nbghostcells !: number of cells in i-direction
INTEGER, PUBLIC :: nbcellsy = jpjglo - 2 - 2*nbghostcells !: number of cells in j-direction
#endif
! local domain size !!! * local computational domain *
INTEGER, PUBLIC :: jpi ! = ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci !: first dimension
INTEGER, PUBLIC :: jpj ! = ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj !: second dimension
INTEGER, PUBLIC :: jpk ! = jpkglo
INTEGER, PUBLIC :: jpim1 ! = jpi-1 !: inner domain indices
INTEGER, PUBLIC :: jpjm1 ! = jpj-1 !: - - -
INTEGER, PUBLIC :: jpkm1 ! = jpk-1 !: - - -
INTEGER, PUBLIC :: jpij ! = jpi*jpj !: jpi x jpj
!!---------------------------------------------------------------------
!! Active tracer parameters
!!---------------------------------------------------------------------
INTEGER, PUBLIC, PARAMETER :: jpts = 2 !: Number of active tracers (=2, i.e. T & S )
INTEGER, PUBLIC, PARAMETER :: jp_tem = 1 !: indice for temperature
INTEGER, PUBLIC, PARAMETER :: jp_sal = 2 !: indice for salinity
INTEGER, PUBLIC, PARAMETER :: jp_dep = 3 !: indice for depth
INTEGER, PUBLIC, PARAMETER :: jp_msk = 4 !: indice for depth
!!----------------------------------------------------------------------
!! Domain decomposition
!!----------------------------------------------------------------------
!! if we dont use massively parallel computer (parameters jpni=jpnj=1) so jpiglo=jpi and jpjglo=jpj
INTEGER, PUBLIC :: jpni !: number of processors following i
INTEGER, PUBLIC :: jpnj !: number of processors following j
INTEGER, PUBLIC :: jpnij !: nb of local domain = nb of processors ( <= jpni x jpnj )
INTEGER, PUBLIC, PARAMETER :: jpr2di = 0 !: number of columns for extra outer halo
INTEGER, PUBLIC, PARAMETER :: jpr2dj = 0 !: number of rows for extra outer halo
INTEGER, PUBLIC, PARAMETER :: jpreci = 1 !: number of columns for overlap
INTEGER, PUBLIC, PARAMETER :: jprecj = 1 !: number of rows for overlap
!!----------------------------------------------------------------------
!! NEMO/OPA 4.0 , NEMO Consortium (2016)
!! $Id: par_oce.F90 7646 2017-02-06 09:25:03Z timgraham $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!======================================================================
END MODULE par_oce
NAMELISTS_AND_FORTRAN_FILES/f_files/sbctide.F90 0 → 100755
View file @ 868014cc
MODULE sbctide
!!======================================================================
!! *** MODULE sbctide ***
!! Initialization of tidal forcing
!!======================================================================
!! History : 9.0 ! 2007 (O. Le Galloudec) Original code
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE dom_oce ! ocean space and time domain
USE phycst ! physical constant
USE daymod ! calandar
USE tideini !
!
USE in_out_manager ! I/O units
USE iom ! xIOs server
USE ioipsl ! NetCDF IPSL library
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
! NB - to access love number
USE bdytides
! END NB
IMPLICIT NONE
PUBLIC
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro !
!!----------------------------------------------------------------------
!! tidal potential
!!----------------------------------------------------------------------
!! sbc_tide :
!! tide_init_potential :
!!----------------------------------------------------------------------
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot, phi_pot
!!----------------------------------------------------------------------
!! NEMO/OPA 3.5 , NEMO Consortium (2013)
!! $Id: sbctide.F90 7646 2017-02-06 09:25:03Z timgraham $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE sbc_tide( kt )
!!----------------------------------------------------------------------
!! *** ROUTINE sbc_tide ***
!!----------------------------------------------------------------------
INTEGER, INTENT( in ) :: kt ! ocean time-step
INTEGER :: jk ! dummy loop index
INTEGER :: nsec_day_orig ! Temporary variable
!!----------------------------------------------------------------------
IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN ! start a new day
!
IF( kt == nit000 ) THEN
ALLOCATE( amp_pot(jpi,jpj,nb_harmo), &
& phi_pot(jpi,jpj,nb_harmo), pot_astro(jpi,jpj) )
ENDIF
!
amp_pot(:,:,:) = 0._wp
phi_pot(:,:,:) = 0._wp
pot_astro(:,:) = 0._wp
!
! If the run does not start from midnight then need to initialise tides
! at the start of the current day (only occurs when kt==nit000)
! Temporarily set nsec_day to beginning of day.
nsec_day_orig = nsec_day
IF ( nsec_day /= NINT(0.5_wp * rdt) ) THEN
kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
nsec_day = NINT(0.5_wp * rdt)
ELSE
kt_tide = kt
ENDIF
CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo )
!
!
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) 'sbc_tide : Update of the components and (re)Init. the potential at kt=', kt
WRITE(numout,*) '~~~~~~~~ '
DO jk = 1, nb_harmo
WRITE(numout,*) Wave(ntide(jk))%cname_tide, utide(jk), ftide(jk), v0tide(jk), omega_tide(jk)
END DO
ENDIF
!
IF( ln_tide_pot ) CALL tide_init_potential
!
! Reset nsec_day
nsec_day = nsec_day_orig
ENDIF
!
END SUBROUTINE sbc_tide
SUBROUTINE tide_init_potential
!!----------------------------------------------------------------------
!! *** ROUTINE tide_init_potential ***
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zcs ! local scalar
!!----------------------------------------------------------------------
DO jk = 1, nb_harmo
!--- NB 11/2017
! love number now provides in tide namelist
zcons = dn_love_number * Wave(ntide(jk))%equitide * ftide(jk)
! ORIGINAL zcons = 0.7_wp * Wave(ntide(jk))%equitide * ftide(jk)
!--- END NB
DO ji = 1, jpi
DO jj = 1, jpj
ztmp1 = amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) )
ztmp2 = -amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,jk) )
zlat = gphit(ji,jj)*rad !! latitude en radian
zlon = glamt(ji,jj)*rad !! longitude en radian
ztmp = v0tide(jk) + utide(jk) + Wave(ntide(jk))%nutide * zlon
! le potentiel est composé des effets des astres:
IF ( Wave(ntide(jk))%nutide == 1 ) THEN ; zcs = zcons * SIN( 2._wp*zlat )
ELSEIF( Wave(ntide(jk))%nutide == 2 ) THEN ; zcs = zcons * COS( zlat )**2
!--- NB 11/2017
! Add tide potential for long period tides
ELSEIF( Wave(ntide(jk))%nutide == 0 ) THEN ; zcs = zcons * (0.5_wp-1.5_wp*SIN(zlat)**2._wp)
!--- END NB
ELSE ; zcs = 0._wp
ENDIF
ztmp1 = ztmp1 + zcs * COS( ztmp )
ztmp2 = ztmp2 - zcs * SIN( ztmp )
zamp = SQRT( ztmp1*ztmp1 + ztmp2*ztmp2 )
amp_pot(ji,jj,jk) = zamp
phi_pot(ji,jj,jk) = ATAN2( -ztmp2 / MAX( 1.e-10_wp , zamp ) , &
& ztmp1 / MAX( 1.e-10_wp, zamp ) )
END DO
END DO
END DO
!
END SUBROUTINE tide_init_potential
!!======================================================================
END MODULE sbctide
NAMELISTS_AND_FORTRAN_FILES/f_files/step.F90 0 → 100755
View file @ 868014cc
MODULE step
!!======================================================================
!! *** MODULE step ***
!! Time-stepping : manager of the ocean, tracer and ice time stepping
!!======================================================================
!! History : OPA ! 1991-03 (G. Madec) Original code
!! - ! 1991-11 (G. Madec)
!! - ! 1992-06 (M. Imbard) add a first output record
!! - ! 1996-04 (G. Madec) introduction of dynspg
!! - ! 1996-04 (M.A. Foujols) introduction of passive tracer
!! 8.0 ! 1997-06 (G. Madec) new architecture of call
!! 8.2 ! 1997-06 (G. Madec, M. Imbard, G. Roullet) free surface
!! - ! 1999-02 (G. Madec, N. Grima) hpg implicit
!! - ! 2000-07 (J-M Molines, M. Imbard) Open Bondary Conditions
!! NEMO 1.0 ! 2002-06 (G. Madec) free form, suppress macro-tasking
!! - ! 2004-08 (C. Talandier) New trends organization
!! - ! 2005-01 (C. Ethe) Add the KPP closure scheme
!! - ! 2005-11 (G. Madec) Reorganisation of tra and dyn calls
!! - ! 2006-01 (L. Debreu, C. Mazauric) Agrif implementation
!! - ! 2006-07 (S. Masson) restart using iom
!! 3.2 ! 2009-02 (G. Madec, R. Benshila) reintroduicing z*-coordinate
!! - ! 2009-06 (S. Masson, G. Madec) TKE restart compatible with key_cpl
!! 3.3 ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
!! - ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
!! 3.4 ! 2011-04 (G. Madec, C. Ethe) Merge of dtatem and dtasal
!! 3.6 ! 2012-07 (J. Simeon, G. Madec. C. Ethe) Online coarsening of outputs
!! 3.6 ! 2014-04 (F. Roquet, G. Madec) New equations of state
!! 3.6 ! 2014-10 (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves
!! 3.7 ! 2014-10 (G. Madec) LDF simplication
!! - ! 2014-12 (G. Madec) remove KPP scheme
!! - ! 2015-11 (J. Chanut) free surface simplification
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! stp : OPA system time-stepping
!!----------------------------------------------------------------------
USE step_oce ! time stepping definition modules
!
USE iom ! xIOs server
IMPLICIT NONE
PRIVATE
PUBLIC stp ! called by nemogcm.F90
!!----------------------------------------------------------------------
!! NEMO/OPA 3.7 , NEMO Consortium (2015)
!! $Id: step.F90 7753 2017-03-03 11:46:59Z mocavero $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
#if defined key_agrif
RECURSIVE SUBROUTINE stp( )
INTEGER :: kstp ! ocean time-step index
#else
SUBROUTINE stp( kstp )
INTEGER, INTENT(in) :: kstp ! ocean time-step index
#endif
!!----------------------------------------------------------------------
!! *** ROUTINE stp ***
!!
!! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.)
!! - Time stepping of LIM (dynamic and thermodynamic eqs.)
!! - Tme stepping of TRC (passive tracer eqs.)
!!
!! ** Method : -1- Update forcings and data
!! -2- Update ocean physics
!! -3- Compute the t and s trends
!! -4- Update t and s
!! -5- Compute the momentum trends
!! -6- Update the horizontal velocity
!! -7- Compute the diagnostics variables (rd,N2, hdiv,w)
!! -8- Outputs and diagnostics
!!----------------------------------------------------------------------
INTEGER :: ji,jj,jk ! dummy loop indice
INTEGER :: indic ! error indicator if < 0
INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt)
!! ---------------------------------------------------------------------
#if defined key_agrif
kstp = nit000 + Agrif_Nb_Step()
IF( lk_agrif_debug ) THEN
IF( Agrif_Root() .and. lwp) WRITE(*,*) '---'
IF(lwp) WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
ENDIF
IF( kstp == nit000 + 1 ) lk_agrif_fstep = .FALSE.
# if defined key_iomput
IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( cxios_context )
# endif
#endif
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! update I/O and calendar
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
indic = 0 ! reset to no error condition
IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
CALL iom_init( cxios_context ) ! for model grid (including passible AGRIF zoom)
IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid
ENDIF
IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init)
CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp
IF( ln_crs ) CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" ) ! tell IOM we are at time step kstp
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( ln_tide ) CALL sbc_tide( kstp ) ! update tide potential
IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
IF( ln_bdy ) CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries
CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice)
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Update stochastic parameters and random T/S fluctuations
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters
IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Ocean physics update (ua, va, tsa used as workspace)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! THERMODYNAMICS
CALL eos_rab( tsb, rab_b ) ! before local thermal/haline expension ratio at T-points
CALL eos_rab( tsn, rab_n ) ! now local thermal/haline expension ratio at T-points
CALL bn2 ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
CALL bn2 ( tsn, rab_n, rn2 ) ! now Brunt-Vaisala frequency
!
! VERTICAL PHYSICS
CALL zdf_bfr( kstp ) ! bottom friction (if quadratic)
! ! Vertical eddy viscosity and diffusivity coefficients
IF( lk_zdfric ) CALL zdf_ric ( kstp ) ! Richardson number dependent Kz
IF( lk_zdftke ) CALL zdf_tke ( kstp ) ! TKE closure scheme for Kz
IF( lk_zdfgls ) CALL zdf_gls ( kstp ) ! GLS closure scheme for Kz
IF( ln_zdfqiao ) CALL zdf_qiao( kstp ) ! Qiao vertical mixing
!
IF( lk_zdfcst ) THEN ! Constant Kz (reset avt, avm[uv] to the background value)
avt (:,:,:) = rn_avt0 * wmask (:,:,:)
avmu(:,:,:) = rn_avm0 * wumask(:,:,:)
avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
ENDIF
IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths
DO jk = 2, nkrnf ; avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk) ; END DO
ENDIF
IF( ln_zdfevd ) CALL zdf_evd( kstp ) ! enhanced vertical eddy diffusivity
IF( lk_zdftmx ) CALL zdf_tmx( kstp ) ! tidal vertical mixing
IF( lk_zdfddm ) CALL zdf_ddm( kstp ) ! double diffusive mixing
CALL zdf_mxl( kstp ) ! mixed layer depth
! write TKE or GLS information in the restart file
IF( lrst_oce .AND. lk_zdftke ) CALL tke_rst( kstp, 'WRITE' )
IF( lrst_oce .AND. lk_zdfgls ) CALL gls_rst( kstp, 'WRITE' )
!
! LATERAL PHYSICS
!
IF( l_ldfslp ) THEN ! slope of lateral mixing
CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density
IF( ln_zps .AND. .NOT. ln_isfcav) &
& CALL zps_hde ( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient
& rhd, gru , grv ) ! of t, s, rd at the last ocean level
IF( ln_zps .AND. ln_isfcav) &
& CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
IF( ln_traldf_triad ) THEN
CALL ldf_slp_triad( kstp ) ! before slope for triad operator
ELSE
CALL ldf_slp ( kstp, rhd, rn2b ) ! before slope for standard operator
ENDIF
ENDIF
! ! eddy diffusivity coeff.
IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp ) ! and/or eiv coeff.
IF( l_ldfdyn_time ) CALL ldf_dyn( kstp ) ! eddy viscosity coeff.
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Ocean dynamics : hdiv, ssh, e3, u, v, w
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor)
IF(.NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors
CALL wzv ( kstp ) ! now cross-level velocity
CALL eos ( tsn, rhd, rhop, gdept_n(:,:,:) ) ! now in situ density for hpg computation
!!jc: fs simplification
!!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636)
!! but ensures reproductible results
!! with previous versions using split-explicit free surface
IF( ln_zps .AND. .NOT. ln_isfcav ) &
& CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient
& rhd, gru , grv ) ! of t, s, rd at the last ocean level
IF( ln_zps .AND. ln_isfcav ) &
& CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
!!jc: fs simplification
ua(:,:,:) = 0._wp ! set dynamics trends to zero
va(:,:,:) = 0._wp
IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) &
& CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment
IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp ) ! bdy damping trends
#if defined key_agrif
IF(.NOT. Agrif_Root()) &
& CALL Agrif_Sponge_dyn ! momentum sponge
#endif
CALL dyn_adv ( kstp ) ! advection (vector or flux form)
CALL dyn_vor ( kstp ) ! vorticity term including Coriolis
CALL dyn_ldf ( kstp ) ! lateral mixing
CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure
CALL dyn_spg ( kstp ) ! surface pressure gradient
! With split-explicit free surface, since now transports have been updated and ssha as well
IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated
CALL div_hor ( kstp ) ! Horizontal divergence (2nd call in time-split case)
IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component)
CALL wzv ( kstp ) ! now cross-level velocity
ENDIF
CALL dyn_bfr ( kstp ) ! bottom friction
CALL dyn_zdf ( kstp ) ! vertical diffusion
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! cool skin
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF ( ln_diurnal ) CALL stp_diurnal( kstp )
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! diagnostics and outputs (ua, va, tsa used as workspace)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats
IF( nn_diacfl == 1 ) CALL dia_cfl( kstp ) ! Courant number diagnostics
IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth)
IF( lk_diadct ) CALL dia_dct( kstp ) ! Transports
CALL dia_ar5( kstp ) ! ar5 diag
IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis
! NB - new harmonic analysis
IF( lk_diaharm_fast ) &
& CALL dia_harm_fast( kstp ) ! Tidal harmonic analysis - restart and faster version
! END NB
CALL dia_wri( kstp ) ! ocean model: outputs
!
IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output
#if defined key_top
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Passive Tracer Model
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL trc_stp ( kstp ) ! time-stepping
#endif
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Active tracers
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
tsa(:,:,:,:) = 0._wp ! set tracer trends to zero
IF( lk_asminc .AND. ln_asmiau .AND. &
& ln_trainc ) CALL tra_asm_inc ( kstp ) ! apply tracer assimilation increment
CALL tra_sbc ( kstp ) ! surface boundary condition
IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr
IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux
IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme
IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends
IF( ln_bdy ) CALL bdy_tra_dmp ( kstp ) ! bdy damping trends
#if defined key_agrif
IF(.NOT. Agrif_Root()) &
& CALL Agrif_Sponge_tra ! tracers sponge
#endif
CALL tra_adv ( kstp ) ! horizontal & vertical advection
CALL tra_ldf ( kstp ) ! lateral mixing
!!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics
!!gm
CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields
IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Set boundary conditions and Swap
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap
!! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.
!! If so:
!! (i) no need to call agrif update at initialization time
!! (ii) no need to update "before" fields
!!
!! Apart from creating new tra_swp/dyn_swp routines, this however:
!! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between
!! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",
!! e.g. a shift of the feedback interface inside child domain.
!! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
!! place.
!!
!!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
CALL tra_nxt ( kstp ) ! finalize (bcs) tracer fields at next time step and swap
CALL dyn_nxt ( kstp ) ! finalize (bcs) velocities at next time step and swap
CALL ssh_swp ( kstp ) ! swap of sea surface height
IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors
!
IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics
!!gm : This does not only concern the dynamics ==>>> add a new title
!!gm2: why ouput restart before AGRIF update?
!!
!!jc: That would be better, but see comment above
!!
IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file
IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters
#if defined key_agrif
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! AGRIF
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL Agrif_Integrate_ChildGrids( stp )
IF( Agrif_NbStepint() == 0 ) THEN ! AGRIF Update
!!jc in fact update is useless at last time step, but do it for global diagnostics
CALL Agrif_Update_Tra() ! Update active tracers
CALL Agrif_Update_Dyn() ! Update momentum
ENDIF
#endif
IF( ln_diaobs ) CALL dia_obs( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Control
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL stp_ctl ( kstp, indic )
IF( indic < 0 ) THEN
CALL ctl_stop( 'step: indic < 0' )
CALL dia_wri_state( 'output.abort', kstp )
ENDIF
!#if defined key_harm_ana
!--- NB Restart for the tidal harmonic analysis
! IF( ln_harm_ana_store ) CALL harm_ana( kstp ) ! Harmonic analysis of tides
!--- END NB -----------------------------------
!# endif
IF( kstp == nit000 ) THEN
CALL iom_close( numror ) ! close input ocean restart file
IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce
IF(lwm.AND.numoni /= -1 ) &
& CALL FLUSH ( numoni ) ! flush output namelist ice (if exist)
ENDIF
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Coupled mode
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!!gm why lk_oasis and not lk_cpl ????
IF( lk_oasis ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges
!
#if defined key_iomput
IF( kstp == nitend .OR. indic < 0 ) THEN
CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
ENDIF
#endif
!
IF( nn_timing == 1 .AND. kstp == nit000 ) CALL timing_reset
!
END SUBROUTINE stp
END MODULE step
NAMELISTS_AND_FORTRAN_FILES/f_files/step_oce.F90 0 → 100755
View file @ 868014cc
MODULE step_oce
!!======================================================================
!! *** MODULE step_oce ***
!! Ocean time-stepping : module used in both initialisation phase and time stepping
!!======================================================================
!! History : 3.3 ! 2010-08 (C. Ethe) Original code - reorganisation of the initial phase
!! 3.7 ! 2014-01 (G. Madec) LDF simplication
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE dom_oce ! ocean space and time domain variables
USE zdf_oce ! ocean vertical physics variables
USE daymod ! calendar (day routine)
USE sbc_oce ! surface boundary condition: ocean
USE sbcmod ! surface boundary condition (sbc routine)
USE sbcrnf ! surface boundary condition: runoff variables
USE sbccpl ! surface boundary condition: coupled formulation (call send at end of step)
USE sbcapr ! surface boundary condition: atmospheric pressure
USE sbctide ! Tide initialisation
USE sbcwave ! Wave intialisation
USE traqsr ! solar radiation penetration (tra_qsr routine)
USE trasbc ! surface boundary condition (tra_sbc routine)
USE trabbc ! bottom boundary condition (tra_bbc routine)
USE trabbl ! bottom boundary layer (tra_bbl routine)
USE tradmp ! internal damping (tra_dmp routine)
USE traadv ! advection scheme control (tra_adv_ctl routine)
USE traldf ! lateral mixing (tra_ldf routine)
USE trazdf ! vertical mixing (tra_zdf routine)
USE tranxt ! time-stepping (tra_nxt routine)
USE tranpc ! non-penetrative convection (tra_npc routine)
USE eosbn2 ! equation of state (eos_bn2 routine)
USE divhor ! horizontal divergence (div_hor routine)
USE dynadv ! advection (dyn_adv routine)
USE dynbfr ! Bottom friction terms (dyn_bfr routine)
USE dynvor ! vorticity term (dyn_vor routine)
USE dynhpg ! hydrostatic pressure grad. (dyn_hpg routine)
USE dynldf ! lateral momentum diffusion (dyn_ldf routine)
USE dynzdf ! vertical diffusion (dyn_zdf routine)
USE dynspg ! surface pressure gradient (dyn_spg routine)
USE dynnxt ! time-stepping (dyn_nxt routine)
USE stopar ! Stochastic parametrization (sto_par routine)
USE stopts
USE bdy_oce , ONLY: ln_bdy
USE bdydta ! open boundary condition data (bdy_dta routine)
USE bdytra ! bdy cond. for tracers (bdy_tra routine)
USE bdydyn3d ! bdy cond. for baroclinic vel. (bdy_dyn3d routine)
USE sshwzv ! vertical velocity and ssh (ssh_nxt routine)
! (ssh_swp routine)
! (wzv routine)
USE domvvl ! variable vertical scale factors (dom_vvl_sf_nxt routine)
! (dom_vvl_sf_swp routine)
USE ldfslp ! iso-neutral slopes (ldf_slp routine)
USE ldfdyn ! lateral eddy viscosity coef. (ldf_dyn routine)
USE ldftra ! lateral eddy diffusive coef. (ldf_tra routine)
USE zdftmx ! tide-induced vertical mixing (zdf_tmx routine)
USE zdfbfr ! bottom friction (zdf_bfr routine)
USE zdftke ! TKE vertical mixing (zdf_tke routine)
USE zdfgls ! GLS vertical mixing (zdf_gls routine)
USE zdfddm ! double diffusion mixing (zdf_ddm routine)
USE zdfevd ! enhanced vertical diffusion (zdf_evd routine)
USE zdfric ! Richardson vertical mixing (zdf_ric routine)
USE zdfmxl ! Mixed-layer depth (zdf_mxl routine)
USE zdfqiao !Qiao module wave induced mixing (zdf_qiao routine)
USE step_diu ! Time stepping for diurnal sst
USE diurnal_bulk ! diurnal SST bulk routines (diurnal_sst_takaya routine)
USE cool_skin ! diurnal cool skin correction (diurnal_sst_coolskin routine)
USE sbc_oce ! surface fluxes
USE zpshde ! partial step: hor. derivative (zps_hde routine)
USE diawri ! Standard run outputs (dia_wri routine)
USE diaptr ! poleward transports (dia_ptr routine)
USE diadct ! sections transports (dia_dct routine)
USE diaar5 ! AR5 diagnosics (dia_ar5 routine)
USE diahth ! thermocline depth (dia_hth routine)
USE diahsb ! heat, salt and volume budgets (dia_hsb routine)
USE diaharm
!--- NB for restart hamonic analysis
USE diaharm_fast ! harmonic analysis of tides (harm_ana routine)
!--- END NB -----------------------------------
USE diacfl
USE flo_oce ! floats variables
USE floats ! floats computation (flo_stp routine)
USE crsfld ! Standard output on coarse grid (crs_fld routine)
USE asminc ! assimilation increments (tra_asm_inc routine)
! (dyn_asm_inc routine)
USE asmbkg
USE stpctl ! time stepping control (stp_ctl routine)
USE restart ! ocean restart (rst_wri routine)
USE prtctl ! Print control (prt_ctl routine)
USE diaobs ! Observation operator
USE in_out_manager ! I/O manager
USE iom !
USE lbclnk
USE timing ! Timing
#if defined key_iomput
USE xios
#endif
#if defined key_agrif
USE agrif_opa_sponge ! Momemtum and tracers sponges
USE agrif_opa_update ! Update (2-way nesting)
#endif
#if defined key_top
USE trcstp ! passive tracer time-stepping (trc_stp routine)
#endif
!!----------------------------------------------------------------------
!! NEMO/OPA 3.7 , NEMO Consortium (2014)
!! $Id: step_oce.F90 7646 2017-02-06 09:25:03Z timgraham $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!======================================================================
END MODULE step_oce
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment