= Changing the CLIO bathymetry and land-sea mask for PI and PMIP4 runs =
[[BR]]

== '''1. PI simulations with an etopo1 bathymetry and land-sea mask (Spring 2019)''' ==

The high resolution topography etopo1 is used to generate a new bathymetry and land-sea mask for CLIO in PI conditions. The default bathymetry is, according to Campin [1997], from etopo5 [1986]. [[BR]]

__Semi-automated pre-processing steps :__ [[BR]]

[[Image(preprocessing_bathychange.2.png)]]

1. The connectivity program writes the chosen variables (mean bathymetry, hypsometry...) into a text file, either on the rotated or regular CLIO grid. It can also produce the connections between ocean basins thanks to the computation of subgrid sills. [[BR]] 
2. The two grids are put together in the bathymetry program, which writes a text file containing the bathymetry with the land-sea mask (0 in every land grid cells). [[BR]] The mask is generated using the hypsometry, a chosen sea-level (-0.5 m for the PI), and a chosen threshold (% of surface of a grid cell above which the cell is defined as ocean - here 40 %). The mask of a few ocean grid cells is manually forced so that all the critical traits stay open. Small isolated seas are closed. [[BR]] The bathymetry is converted into vertical levels using the bath.om file, which contains the old depth as a vertical level (from 1 to 20, 20 being the deepest cell) and the definition of these irregular vertical levels (center and height). The old vertical level is forced for a few problematic grid cells in order to get realistic salinity values in the Mediterranean Sea and Hudson Bay. The vertical level 1 is avoided (either forced to 0 or 2), because the model cannot deal with these very shallow grid cells. As the model also cannot deal with 'isolated' oceanic grid cells for which the last vertical level is isolated (e.g. deep grid cells with shallower neighbours), a process similar to a smoothing filter is applied. [[BR]] 
3. The pre-poids program loads the new bathymetry (bath.txt file) and generate two files necessary to the next step. It contains the defgrid.f routine, slightly modified to deal the newly generated 'isolated' cells at the edges (longitudes 2 and 121). [[BR]]
4. The poids program uses the clio.msk and clio.srf files to compute the interpolation points between the CLIO and the T21 grid (mozaic.w), as well as the fraction of ocean seen by each T21 grid cells (fractoc.dat, aka fort.8). When running this program, kamax = 14 and komax = 17 should be printed. [[BR]]

__How to run a PI simulation with the new bathymetry and land-sea mask :__ [[BR]]

Run with -F -I (-P), using a scenario (pi-etopo1 - commited but outdated) containing :
- inputdata : fractoc.dat, clio/bath.txt, clio/mozaic.w, clio/bath.om (definition of vertical levels still needed)
- parameters : run.param, with a BERING_VALUE tuned to get a realistic Bering transport (around 0.8 Sv)
- sources : defgrid.f (with the right bath.txt name and modification for isolated cells at the edges)

__Potential sources of error when running :__ [[BR]]
- "no ocean basin defined in x,y" : change the ocbas.dat file in ecbilt_clio/parameters or include a modified one in scenario/inputdata.

- "error in runoff budget" : set bricole_Caspienne=.false. in ludus-code/lbm/sources/routageEAU.f90.

- "Albedo out of range" : conflict with alb.dat, comment the LGMSWITCH.EQ.1 condition in landmodel0.f for now.

- "rain larger than maximum set by rainmax" : you're in trouble, variables are diverging (e.g. due to a non-zero sea ice fraction in a few tropical grid cells - which were formerly land but are now ocean ?). Check if the number of 'isolated' grid cells has fallen to 0 in step 2.

- NaN : try a restart.

== '''2. LGM simulations with a PMIP4 bathymetry and land-sea mask (Spring 2019)''' ==

'''2.1. With Peltier ICE-6G-C (LGM4_P1.1)''' [[BR]]

__Semi-automated pre-processing :__ [[BR]]

The same four pre-processing steps are applied with the following LGM topography, computed with a regridding of the Peltier anomaly on the etopo1 grid (step 0) : [[BR]]
{{{
LGM bathymetry = PI (etopo1) + LGM Peltier (ICE-6G-C 21kyr) - PI Peltier (ICE-6G-C 0kyr)
}}}

''Note '' : For the LGM, the chosen sea-level (in step 2) can either be consistent with Lambeck data (-133.9 m for the LGM) or with the ice sheets reconstructions of Peltier or Tarasov. Because of this lower sea-level, the manual points where the land-sea mask is forced have to be redefined. Some stay the same (Gibraltar), others are not necessary anymore (Hudson Bay and Japan Sea outlets) while a few new critical points appear (Fram Strait, Golf of Mexico outlet). [[BR]]

__How to run a LGM PMIP4 simulation with the new bathymetry and land-sea mask :__ [[BR]]

Run with -F -I -P, using a scenario (lgm-pmip4-peltier, commited but outdated) also containing : [[BR]]
- inputdata : fractoc.dat, clio/bath.txt, clio/mozaic.w, clio/bath.om
- parameters : run.param, with BERING_VALUE set to 0
- sources : defgrid.f
As well as the following files :
- sources : choixcomposantes.h with LGMSWITCH 1, F_PALAEO 3 and ISM 3
- inputdata : berg.dat and icemask.dat (necessary for the initialization but regenerated thanks to ISM 3)
- sources : transfer_ecb_grisli_mod.f90 (to write less outputs) and input_GRISLI2L.f90 -> input_GRISLI2L.f90.peltier26k [[BR]]
''Note '' : With F_PALAEO 3, don't forget to set palaeo_year = 21000 in ludus-code/ecbilt/parameters/palaeo_transient_parameters.dat [[BR]]

'''IMPORTANT :''' Set the symbolic link in ludus-code/ecbilt/inputdata/ : Gano_40k-0k_hemin40.nc -> ICE6G_26k-0k_hemin40.nc [[BR]]


'''2.2. With Tarasov GLAC-1D (LGM4_T1.1)''' [[BR]]

The same pre-processing steps are applied to generate the bath.txt, fractoc.dat and mozaic.w files, including step 0 :

{{{
LGM bathymetry = PI (etopo1) + LGM Tarasov (GLAC-1D 21kyr) - PI Tarasov (GLAC-1D 0kyr)
}}}

''Note '' : The manual points defined to force the land-sea mask in critical straits are very similar to the ones defined for the LGM4_P1.1. [[BR]]

Run with -F -I -P, using a scenario (lgm-pmip4-tarasov, commited but outdated) containing these three files, along with the other files mentioned in the previous section. The same files are contained in sources, except for input_GRISLI2L.f90 (-> input_GRISLI2L.f90.tarasov26k).

'''IMPORTANT :''' Set the symbolic link in ludus-code/ecbilt/inputdata/ : Gano_40k-0k_hemin40.nc -> GLAC1D_26k-0k_hemin40.nc
 [[BR]]

== '''3. Simulations with all the PMIP4 boundary conditions (Fall 2020)''' ==

Following the implementation of the icemask and elevation as netcdf files, we have properly generated a berg.nc and icemask.nc corresponding to either the ICE-6G-C or the GLAC-1D ice-sheet reconstruction, instead of using old berg.dat and icemask.dat files. As a result, the PMIP4 elevation associated with the up-to-date reconstructions is implemented not only in the Northern Hemisphere (old trick using ISM 3) but also in Antarctica. [[BR]] [[BR]]
To do that, the variables extracted from Peltier's or Tarasov's topography are : [[BR]]
1. interpolated on the ECBilt grid using cdo remapbil [[BR]]
2. (for berg.nc) set to 0 if the topography is negative (oceans) with : 
{{{
ncap2 -v -O -s agg1=agg1; where (agg1<0) agg1=0 orogT21.nc orogT21_agg1corr.nc
}}}
3. (for berg.nc) attached to the agg2 variable (which we don't know how to regenerate properly) :
{{{
ncks -A orogT21_agg1corr.nc agg2.nc
}}}
2. (for icemask.nc) divided by 100 to get values between 0 and 1 (only for ICE-6G-C) [[BR]]
3. (for icemask.nc) split to either 0 or 1 (as the interpolation creates values in-between), using a threshold value of 20% to ensure that the surface area covered by ice is as close to the one from reconstructions as possible : [[BR]]
{{{
		if icemask[i,j]>=0.20:
			icemask[i,j]=1
		else:
			icemask[i,j]=0
}}}
4. additioned/subtracted to get the following anomaly : [[BR]]

{{{
LGM topography = PI (etopo1) + LGM Peltier (ICE-6G-C 21kyr) - PI Peltier (ICE-6G-C 0kyr)
}}}

Run with -F -I -P, using a scenario such as lgm-pmip4-fullICE6GC-r1344 (revision 1355) which contains : [[BR]]
- inputdata : fractoc.dat, clio/bath.txt, clio/mozaic.w, but also icemask.nc and berg.nc
- parameters : run.param, with BERING_VALUE set to 0
- sources : BC_switches with NC_IMSK 1 and NC_BERG 1, and choixcomposantes.h with LGMSWITCH 1, F_PALAEO 2, ISM 0, and BATHY 1

''Note '' : No need for the symbolic link in ludus-code/ecbilt/inputdata/ anymore. But don't forget to set palaeo_year = 21000 in ludus-code/ecbilt/parameters/palaeo_transient_parameters.dat  [[BR]]

== '''4. PMIP4 Transient simulations (see Nathaëlle's page)''' ==

To accomodate restarts with a modified bathymetry (and in the long run an interactive bathymetry for transient simulations), the model has to generate realistic values of the tracers content of new oceanic grid cells. A new module has been commited (clio/sources/update_clio_bathy_tools.f90). The initialization of all the restart variables in new ocean grid cells (in redrunb.f) is now done by either setting 0 or averaging the values in neighbouring oceanic grid cells.