Updated 2020-02-20
The following sketch displays our hierarchical experiments trying to understand the origin of large-scale circulations.
A series of articles archived the procedures to obtain the implementation of the experiments:
Now, finally, we chose the red pill, try to go deep into the model configuration/source, to change the land/sea mask and implement our pure-aqua and surf-aqua experiments.
The implementation of fully coupled aqua-planet experiment is not so straightforward since there is only official support for the prescribed SST aqua-planet. Therefore, as deep-time paleo simulation need to change the land/sea distribution, we seek to use the paleo toolkit to achieve our goal. See reference.
We first begin with the ocean model, POP2. In the guide for deep-time simulation:
The deep time paleo modeler is responsible for creating the binary ocean bottom topography file, the horizontal and vertical ocean grid files, a region mask file, and the ocean initial conditions. In addition, several input template files may need to be changed. For example, if the land/ocean mask is changed, the input template file containing new indices for diagnostic transport calculations will need to be changed. The namelist input file (pop2_in and user_nl_pop2) will need to be changed as well. This section will describe methods that can be used to generate these new files.
Okay, there are four files need to be created.
In pop2_in
, we see the namelist variable define oceanic bottom topography:
topography_file = '/users/yangsong3/CESM/input/ocn/pop/gx1v6/grid/topography_20090204.ieeei4'
This is an IEEE standard 4 byte integer binary file. Since the paleo toolkit for changing the bathymetry is quite complex, we tried to use basic ncl script to change the value.
Results from the file:
(0) min=0 max=60
60 means the deepest ocean bottom lies on the 60th model layer. A list of the ocean depth is attached below:
1: 10m; 2: 20m; 3: 30m; 4: 40m; 5: 50m; 6: 60m; 7: 70m; 8: 80m; 9: 90m; 10: 100m; 11: 110m; 12: 120m; 13: 130m; 14: 140m; 15: 150m; 16: 160m; 17: 170.1968m; 18: 180.7613m; 19: 191.8212m; 20: 203.4993m; 21: 215.9234m; 22: 229.2331m; 23: 243.5845m; 24: 259.1558m; 25: 276.1526m; 26: 294.8147m; 27: 315.4237m; 28: 338.3123m; 29: 363.8747m; 30: 392.5805m; 31: 424.9889m; 32: 461.7666m; 33: 503.7069m; 34: 551.7491m; 35: 606.9967m; 36: 670.7286m; 37: 744.398m; 38: 829.607m; 39: 928.0435m; 40: 1041.368m; 41: 1171.04m; 42: 1318.094m; 43: 1482.901m; 44: 1664.992m; 45: 1863.014m; 46: 2074.874m; 47: 2298.039m; 48: 2529.904m; 49: 2768.098m; 50: 3010.671m; 51: 3256.138m; 52: 3503.449m; 53: 3751.892m; 54: 4001.012m; 55: 4250.525m; 56: 4500.261m; 57: 4750.12m; 58: 5000.047m; 59: 5250.009m; 60: 5499.991m;
Here we set all to Layer 25 49 in PURE_AQUA
(2.768km depth), and all land (0) to 1 (10m) in SURF_AQUA
model.
However, the NCL had problem when writing the new file (See Appendix01 for details and solution). We change to F90 to realize our goal. For example, below is the code to generate SURF_AQUA
bathymetry data:
program change_bathymetry
implicit none
integer, parameter :: fin = 100, fout= 101
integer, parameter :: nlen = 122880
integer :: ii
integer :: status = 0
integer :: bath(nlen)=0
character(len=256) :: filename="/users/yangsong3/CESM/input/ocn/pop/gx1v6/grid/topography_20090204.ieeei4"
character(len=256) :: outname="/users/yangsong3/CESM/input/ocn/pop/gx1v6/grid/topography_surf_aqua_20180607.ieeei4"
open(unit=fin, file=filename, form="unformatted", access="direct", recl=nlen, convert="big_endian")
read(fin, rec=1) bath
close(fin)
do ii = 1, nlen
if (bath(ii)<1) then
bath(ii)=1
end if
end do
open(unit=fout, file=outname, form="unformatted", access="direct", action="write", recl=nlen, convert="big_endian")
write(fout, rec=1) bath
close(fout)
end program
Using the above program, we rewrote two files with SURF_AQUA
and PURE_AQUA
settings.
We first check the horizontal grid system.
It can be seen that over the Northern Polar region, the grid system shifts the north pole onto the Greenland Island, and both NH and SH grid poles are located over the landmass.
According to the Paleo FAQ:
The ocean model requires that grid poles be placed over land. Numerically no computation can be done at the convergence point of all longitudes at the grid pole. The ocean model solves this problem by shifting the grid pole away from the geographic pole and placing it over a land mass. (Atmospheric models solve this problem by using numerical filters). Therefore if there is no land at the geographic pole, the numerical pole must be shifted over land elsewhere.
To make it simple, we first tried to maintain the default Greenland replaced grid pole/ SH pole with no land mass. If that does not work, we will try to maintain the Greenland and Antarctica landmass, which will not influence our main results.
(The upper test does not work! We have to maintain the polar land inside the ocean grid “black hole”! Otherwise, at the first step of integration, the atmosphere will complain it got a zero longwave radiation (temperature=0K) from the black hole.).
We check the POP default basin mask file, with the min=-14 and max=11, the masking figure can be found on the NCL popmask page.
And the gx1v6_region_idx
is avaliable in the run folder:
IDX | OCEAN NAME |
---|---|
1 | ‘Southern Ocean ‘ |
2 | ‘Pacific Ocean ‘ |
3 | ‘Indian Ocean ‘ |
4 | ‘Persian Gulf ‘ |
-5 | ‘Red Sea ‘ |
6 | ‘Atlantic Ocean ‘ |
7 | ‘Mediterranean Sea ‘ |
8 | ‘Labrador Sea ‘ |
9 | ‘GIN Sea ‘ |
10 | ‘Arctic Ocean ‘ |
11 | ‘Hudson Bay ‘ |
-12 | ‘Baltic Sea ‘ |
-13 | ‘Black Sea ‘ |
-14 | ‘Caspian Sea ‘ |
Minus-zero Seas are for marginal seas.
We then follow the guide for Deep Time simulation to change the idx:
Deep Time: Since deep time simulations have ocean configurations that differ profoundly from current day, deep time modelers will need to define their own ocean regions. Typically we recommend that they define two simple regions, Northern and Southern Hemisphere.
We change gx1v6_region_idx
to:
IDX | OCEAN NAME |
---|---|
1 | ‘Southern Ocean ‘ |
2 | ‘Northern Ocean ‘ |
region_mask_gx1v6
is modified accordingly.
The final work is to change the ocean initial condition. For PURE_AQUA
and SURF_AQUA
, the zonal-mean type should fit well.
ZONAL-MEAN
Initialize with a global, zonally averaged temperature/salinity distribution. POP allows a zonally averaged temperature/salinity file to be used for initialization. This file is binary and the format can be found in the CCSM3/POP source code, subroutine initial.F
Simply change the pop2_in
:
init_ts_option = "zonal-mean"
However, it seems that this is a legacy configuration in CCSM3/POP. In POP2, I cannot find the zonal-mean option in the source code, so I conducted a simple test to prove this. Unfortunately, the model does not recognize the settings. Therefore, I have to change the default initial condition to fit the aqua-planet run.
The default initial file for temperature and salinity is ts_PHC2_jan_ic_gx1v6_20090205.ieeer8
. Thus, in the SURF_AQUA
experiment, we use the poisson grid fill to fill the grid on the ground, and in the PURE_AQUA
experiment, we fill all grid with zonal mean and mask out layer deeper than L26.
We should also use the recommended namelist settings for Deep Time Simulations to turn off the processes not so related to the climate physics:
overflows_on = .false.
overflows_interactive = .false.
lhoriz_varying_bckgrnd = .false.
ldiag_velocity = .false.
ltidal_mixing = .false.
bckgrnd_vdc1 = 0.524
bckgrnd_vdc2 = 0.313
sw_absorption_type = 'jerlov'
jerlov_water_type = 3
chl_option = 'file'
chl_filename = 'unknown-chl'
chl_file_fmt = 'bin'
All CESM Land/Sea mask info comes from the ocean data. In the paleo resources, we follow the instrudctions to make our mapping and domain data.
Notes:
${CESM_INPUT}/share/scripgrids/
or download heregx1v6_090205.nc
grid_imask variable ${CESMROOT}/tools/mapping/gen_mapping_files/gen_cesm_maps.sh -fatm /users/yangsong3/CESM/input/share/scripgrids/fv1.9x2.5_090205.nc -natm fv19_25 -focn /users/yangsong3/CESM/input/share/scripgrids/gx1v6_aqua_polar_180614.nc -nocn gx1PT --nogridcheck
/users/yangsong3/CESM/cesm1_2_2/tools/mapping/gen_domain_files
to build gen_domain
executable. Note to change the following line in Macro if met MPI/NC variable not found problem:
SFC:= mpif90
./gen_domain -m ../gen_mapping_files/map_gx1PT_TO_fv19_25_aave.180614.nc -o gx1PT -l fv19_25
Similar as Coupler generated mapping files, there is a land-sea mask file /users/yangsong3/CESM/input/atm/cam/ocnfrac/domain.camocn.1.9x2.5_gx1v6_090403.nc
, using similar technique, change all mask
and frac
to 1 to present aqua-planet.
Force the model generate sea ice by itself:
Ice initial conditions files are not required for deep time paleoclimate cases. We recommend initializing the ice model with a ‘no ice’ initial state and allowing the model to simulate an ice state. Set ‘no ice’ in the ice model namelist (user_nl_cice).
ice_ic = 'none'
Please see Appendix02 for detailed namelist changes.
Land model never works in the aqua-planet, so we just use the cpl-generated domain file to drive the model. However, in NO_TOPO
series experiments, we reset the surface data to bare ground.
We first found that when we use the NCL to write a direct binary file, the file records very large values near the top of the variable type (e.g. 21474836xx for int, and 3e317 for double). Therefore, we tried to use fortran to deal with the problem. Finally, we found this is actually we set the file option as Big_Endian in NCL only for read process. Solution is simple, also set for write process:
setfileoption("bin","WriteByteOrder","BigEndian")
Note to use WriteByteOrder
.
<!--"atm domain file (char) " -->
<entry id="ATM_DOMAIN_FILE" value="domain.lnd.fv1.9x2.5_gx1v6.aqua.polar.180614.nc" />
<!--"lnd domain file (char) " -->
<entry id="LND_DOMAIN_FILE" value="domain.lnd.fv1.9x2.5_gx1v6.aqua.polar.180614.nc" />
<!--"ice domain file (char) " -->
<entry id="ICE_DOMAIN_FILE" value="domain.ocn.gx1v6.aqua.polar.180614.nc" />
<!--"ocn domain file (char) " -->
<entry id="OCN_DOMAIN_FILE" value="domain.ocn.gx1v6.aqua.polar.180614.nc" />
<!--"atm2ocn flux mapping file (char) " -->
<entry id="ATM2OCN_FMAPNAME" value="cpl/gridmaps/fv1.9x2.5/map_fv19_25_TO_gx1v6_aave.aqua.polar.180614.nc" />
<!--"atm2ocn state mapping file (char) " -->
<entry id="ATM2OCN_SMAPNAME" value="cpl/gridmaps/fv1.9x2.5/map_fv19_25_TO_gx1v6_blin.aqua.polar.180614.nc" />
<!--"atm2ocn vector mapping file (char) " -->
<entry id="ATM2OCN_VMAPNAME" value="cpl/gridmaps/fv1.9x2.5/map_fv19_25_TO_gx1v6_patc.aqua.polar.180614.nc" />
<!--"ocn2atm flux mapping file (char) " -->
<entry id="OCN2ATM_FMAPNAME" value="cpl/gridmaps/gx1v6/map_gx1v6_TO_fv19_25_aave.aqua.polar.180614.nc" />
<!--"ocn2atm state mapping file (char) " -->
<entry id="OCN2ATM_SMAPNAME" value="cpl/gridmaps/gx1v6/map_gx1v6_TO_fv19_25_aave.aqua.polar.180614.nc" />
init_ts_file = '/users/yangsong3/CESM/input/ocn/pop/gx1v6/ic/ts_surf_aqua_polar_PHC2_jan_ic_gx1v6_20180613.ieeer8'
region_mask_file = '/users/yangsong3/CESM/input/ocn/pop/gx1v6/grid/region_mask_aqua_polar_20180612.ieeei4'
region_info_file = '/users/yangsong3/CESM/cesm1_2_2/models/ocn/pop2/input_templates/gx1v6_aqua_region_ids'
topography_file = '/users/yangsong3/CESM/input/ocn/pop/gx1v6/grid/topography_surf_aqua_polar_20180607.ieeei4'
n_transport_reg = 1
diag_transp_freq_opt='never'
diag_transport_file='/users/yangsong3/CESM/cesm1_2_2/models/ocn/pop2/input_templates/gx1v6_transport_contents_aqua'
!-- dt_count sets the number of times the ocean is coupled per NCPL_BASE_PERIOD.
dt_count=30
overflows_on = .false.
overflows_interactive = .false.
lhoriz_varying_bckgrnd = .false.
ldiag_velocity = .false.
ltidal_mixing = .false.
!--The new background vertical velocity parameters were determined by Gokhan and Christine to increase Kappa-v in the deeper ocean in absence of tidal mixing. The tidal mixing normally would do this.
bckgrnd_vdc1 = 0.524
bckgrnd_vdc2 = 0.313
sw_absorption_type = 'jerlov'
jerlov_water_type = 3
chl_option = 'file'
chl_filename = 'unknown-chl'
chl_file_fmt = 'bin'
bnd_topo = '/users/yangsong3/CESM/input/atm/cam/topo/USGS-gtopo30_1.9x2.5_pure.nc'
ncdata='/users/yangsong3/L_Zealot/B/B20f19-surf-aqua/indata/B20f19-topo.cam.i.0010-01-01-00000.nc'
!dtime=900
finidat = '/users/yangsong3/L_Zealot/B/B20f19-surf-aqua/indata/B20f19-surf-aqua.clm2.r.0001-01-06-00000.nc'
fsurdat = '/users/yangsong3/CESM/input/lnd/clm2/surfdata/surfdata_1.9x2.5_bare_ground_simyr2000_c091005.nc'
urban_hac='OFF'
kmt_file = '/users/yangsong3/CESM/input/ocn/pop/gx1v6/grid/topography_surf_aqua_polar_20180607.ieeei4'
ice_ic = 'none'
xndt_dyn = 2
finidat_rtm=''
POP aborting… (init_moc_ts_transport_arrays) no transport regions have been detected – che ck namelist transports_nml
It seems the ocean transport diag met some problems.
The transport_contents input_template is used for diagnostic purposes only. In this file, the modeler can specify ocean locations (straits for example) for ocean transport computations that will be output in the ocean log files.
I first tried to create my own transport_contents
file, and we met end-of-file read crashes. It is clear that the number of diag locations is hard-coded in the model source code. So I then tried to change the namelist variable:
diag_transp_freq_opt = 'never'
The same error message as the beginning… Then I found the first line in the transport_contents
file denotes the number of diag places. So I change it to zero. Another time, not work.
I then found another namelist variable:
n_transport_reg=1
Originally, 2, and I change to 1. New error occurs!
(seq_map_avNormAvF) ERROR: lsize_i ne lsize_f 0 245209 MCT::m_MatAttrVectMul::sMatAvMult_SMPlus_:: FATAL ERROR–parallelization strategy name not supported. 000.MCT(MPEU)::die.: from MCT::m_MatAttrVectMul::sMatAvMult_SMPlus_()
Oops! Big problem in the decomposition part to parallel computing.
After long-time test and turn on the detailed debug info, we got another error indicating
forrtl: error (73): floating divide by zero
Then I found the error occurs in the radiation module. This is because when the atm model compute radiation, it is faced with absolute zero inside polar region where the ocean grid cannot cover.
Thus, we have to maintain polar land inside the ocean grid “black hole” to make the model run.
In the PURE_AQUA
experiment, the model crashes after several months’ integration. I have tried to change the initial condition, increase the timesteps, and modify the oceanic topography to add a smooth change from polar land to 2.7km depth ocean. All the trials failed. Finally, when I remain the sharp change from polar land to 300m depth water, and remain the oceanic topography deeper than 300m, the model runs smoothly. I have no idea why this happens but it seems to be a very wiered problem. Anyway, 300m-depth mixed layer is enough for our scientific question.
Ocean
Coupler
Updated 2020-02-20
A new NCC paper provided a technology to rectify the fully-coupled model systematic bias.
The procedures are following:
FREE
) after spin-up, of course, full of systematic biases.SSTR
) toward observational SST.
SSTR-FREE
, make it as a constant forcing and add the forcing to free run, namely HEATR_CTRL
.HEATR_SEN
).I will then show the procedures to conduct the HEATR_CTRL
run.
First we lock the definition part for the cam variables to hold data input to & output from the coupler. In $CESMROOT/models/atm/cam/src/control
:
!
! Public data types
!
public cam_out_t ! Data from atmosphere
public cam_in_t ! Merged surface data
!---------------------------------------------------------------------------
! This is the data that is sent from the atmosphere to the surface models
!---------------------------------------------------------------------------
type cam_out_t
integer :: lchnk ! chunk index
integer :: ncol ! number of columns in chunk
real(r8) :: tbot(pcols) ! bot level temperature
real(r8) :: zbot(pcols) ! bot level height above surface
real(r8) :: ubot(pcols) ! bot level u wind
real(r8) :: vbot(pcols) ! bot level v wind
real(r8) :: qbot(pcols,pcnst) ! bot level specific humidity
real(r8) :: pbot(pcols) ! bot level pressure
real(r8) :: rho(pcols) ! bot level density
real(r8) :: netsw(pcols) !
real(r8) :: flwds(pcols) !
real(r8) :: precsc(pcols) !
real(r8) :: precsl(pcols) !
real(r8) :: precc(pcols) !
real(r8) :: precl(pcols) !
real(r8) :: soll(pcols) !
real(r8) :: sols(pcols) !
real(r8) :: solld(pcols) !
real(r8) :: solsd(pcols) !
real(r8) :: thbot(pcols) !
real(r8) :: co2prog(pcols) ! prognostic co2
real(r8) :: co2diag(pcols) ! diagnostic co2
real(r8) :: psl(pcols)
real(r8) :: bcphiwet(pcols) ! wet deposition of hydrophilic black carbon
real(r8) :: bcphidry(pcols) ! dry deposition of hydrophilic black carbon
real(r8) :: bcphodry(pcols) ! dry deposition of hydrophobic black carbon
real(r8) :: ocphiwet(pcols) ! wet deposition of hydrophilic organic carbon
real(r8) :: ocphidry(pcols) ! dry deposition of hydrophilic organic carbon
real(r8) :: ocphodry(pcols) ! dry deposition of hydrophobic organic carbon
real(r8) :: dstwet1(pcols) ! wet deposition of dust (bin1)
real(r8) :: dstdry1(pcols) ! dry deposition of dust (bin1)
real(r8) :: dstwet2(pcols) ! wet deposition of dust (bin2)
real(r8) :: dstdry2(pcols) ! dry deposition of dust (bin2)
real(r8) :: dstwet3(pcols) ! wet deposition of dust (bin3)
real(r8) :: dstdry3(pcols) ! dry deposition of dust (bin3)
real(r8) :: dstwet4(pcols) ! wet deposition of dust (bin4)
real(r8) :: dstdry4(pcols) ! dry deposition of dust (bin4)
end type cam_out_t
!---------------------------------------------------------------------------
! This is the merged state of sea-ice, land and ocean surface parameterizations
!---------------------------------------------------------------------------
type cam_in_t
integer :: lchnk ! chunk index
integer :: ncol ! number of active columns
real(r8) :: asdir(pcols) ! albedo: shortwave, direct
real(r8) :: asdif(pcols) ! albedo: shortwave, diffuse
real(r8) :: aldir(pcols) ! albedo: longwave, direct
real(r8) :: aldif(pcols) ! albedo: longwave, diffuse
real(r8) :: lwup(pcols) ! longwave up radiative flux
real(r8) :: lhf(pcols) ! latent heat flux
real(r8) :: shf(pcols) ! sensible heat flux
real(r8) :: wsx(pcols) ! surface u-stress (N)
real(r8) :: wsy(pcols) ! surface v-stress (N)
real(r8) :: tref(pcols) ! ref height surface air temp
real(r8) :: qref(pcols) ! ref height specific humidity
real(r8) :: u10(pcols) ! 10m wind speed
real(r8) :: ts(pcols) ! merged surface temp
real(r8) :: sst(pcols) ! sea surface temp
real(r8) :: snowhland(pcols) ! snow depth (liquid water equivalent) over land
real(r8) :: snowhice(pcols) ! snow depth over ice
real(r8) :: fco2_lnd(pcols) ! co2 flux from lnd
real(r8) :: fco2_ocn(pcols) ! co2 flux from ocn
real(r8) :: fdms(pcols) ! dms flux
real(r8) :: landfrac(pcols) ! land area fraction
real(r8) :: icefrac(pcols) ! sea-ice areal fraction
real(r8) :: ocnfrac(pcols) ! ocean areal fraction
real(r8), pointer, dimension(:) :: ram1 !aerodynamical resistance (s/m) (pcols)
real(r8), pointer, dimension(:) :: fv !friction velocity (m/s) (pcols)
real(r8), pointer, dimension(:) :: soilw !volumetric soil water (m3/m3)
real(r8) :: cflx(pcols,pcnst) ! constituent flux (evap)
real(r8) :: ustar(pcols) ! atm/ocn saved version of ustar
real(r8) :: re(pcols) ! atm/ocn saved version of re
real(r8) :: ssq(pcols) ! atm/ocn saved version of ssq
real(r8), pointer, dimension(:,:) :: depvel ! deposition velocities
end type cam_in_t
Here we note that SW and LW are from coupler, and LH and SH are from coupler. Then we will grep the relationship between these things and the variables in the output file of CAM model.
Name in code | Name in output |
---|---|
cam_out%netsw | FSNS |
cam_out%flwds | FLDS |
cam_in%lhf | LHFLX |
cam_in%shf | SHFLX |
Next, in physpkg.F90
(still this guy!), we found that lhf is calculated by cam_in%cflx (moisture flux from the surface). As the moisture flux is a constituent flux, I decide to remain the latent heat at first.
Therefore, the question is to add a constant forcing to netsw, flwds, and shflx.
Luckily, we have tried to change the shflx in previous trial.
For netsw and flwds, physpkg.F90
call the radiation, and in radiation_tend module, the netsw is set. Note that in radiation.F90
:
real(r8), parameter :: cgs2mks = 1.e-3_r8
flds(i) = cam_out%flwds(i)*cgs2mks
cam_out%flwds(i) = cam_out%flwds(i)*cgs2mks
I don’t know why they organize the data structure like this. But netsw and flwds both can be handled similarly after physpkg call the radiation module.
The modified code files have been uploaded to github.
Updated 2018-11-01
开始测试Aqua-Planet试验,直接按照aqua的设置run没有任何问题,但是发现SST没有annual cycle,询问师姐:
你是用的cesm的default的 aqua compset吗 如果是的话,那么好像是不需要额外读入sst文件的,而是直接在code里指定需要哪种sst profile就可以,具体是通过修改ocn_comp.F90 中 sst_option,这里不同的sst_option代表不同的profile;关于topo,好像你是需要在user_nl_cam里指定一个新的topo文件的(把高度啥的设为0)
这么说还需要自己写一个topo为0的文件。
修改轨道参数,固定黄赤交角为0度:
rb_mode='fixed_parameters'
orb_obliq=0.0
orb_eccen=0.0167239
orb_mvelp=102.035
Updated 2018-05-14