Please check my GitHub Repo for the most recent code implementation in CESM1.2.2.
在之前MM高原试验基础上，应黄老师要求增加北大西洋SST restoring的试验。因通量计算在海洋模块进行，所以应该在海洋模块计算之前修改并update。
查到论坛帖子，以及HC师兄提供的code.
看帖子似乎有namelist变量控制。找到一个NCAR的github主页内容,不过是针对只转海洋模式的情况。
经查HC师兄的code，发现思路比较简单，就是在海洋模块每一次积分step后，用人工读入的文件做restoring。不过貌似还有可以优化的地方，读文件模块应该可以放到其他module，避免模式每一积分步都在读写。
测试编译一次通过。Nice.关键是拿到对应的文件，并进一步优化一下程序流程。
POP的history文件比较诡异，不仅有.h文件，还有.h.nday1文件。
在faq找到了关于nday1的说明
By default, the daily averaged fields are bundled into monthly timeseries files. These files are relatively small, because they contain only a handful of 2D fields. Out of the box, these fields are: Surface Potential Temperature (SST), SST2, Mixed-Layer Depth, and Maximum Mixed-Layer Depth.
这个文件在之前run ctrl的时候并没有存下，那么接下来的问题就是月平均文件是否足以提供full layer海温，包括SST呢？
瞄了一下h文件，发现第一层代表10m深的水，那么应该是没问题的了，直接restore 10m水深。
Updated 2018-01-28
继续看，存储位温的变量是TRACER，但这个变量比较诡异，是一个6维的数据，0 1 2 容易理解，是经纬度，和层次数，第1层为10m水深。3不知道是什么，4是previous/current时间帧，最后一维是iblock用于MPI分片。 问题就是神奇的第3维。
检查了下程序，发现第3维的循环控制变量bound为nt，step_mod.F90内有以下语句：
SBUFF_SUM(:,:,iblock,index_o2x_So_t ) = &
SBUFF_SUM(:,:,iblock,index_o2x_So_t ) + delt* &
TRACER(:,:,1,1,curtime,iblock)
SBUFF_SUM(:,:,iblock,index_o2x_So_s ) = &
SBUFF_SUM(:,:,iblock,index_o2x_So_s ) + delt* &
TRACER(:,:,1,2,curtime,iblock)
进一步查找，发现prognostic.F90有以下定义
type (tracer_field), dimension(nt) :: &
tracer_d ! descriptors for each tracer
!......
tracer_d(1)%short_name = 'TEMP'
tracer_d(1)%long_name = 'Potential temperature'
tracer_d(1)%units = 'degC'
tracer_d(1)%tend_units = 'degC/s'
tracer_d(1)%flux_units = 'degC cm/s'
tracer_d(1)%scale_factor = 1.0_rtavg
tracer_d(2)%short_name = 'SALT'
tracer_d(2)%long_name = 'Salinity'
tracer_d(2)%units = 'msu (g/g)'
tracer_d(2)%tend_units = 'gram/kilogram/s'
tracer_d(2)%flux_units = 'gram/kilogram cm/s'
tracer_d(2)%scale_factor = 1000.0_rtavg
明白啦，1是位温，2是盐度。所以修改思路就很明确了：
- 读入SST文件和weight文件
- 在每一个step后对SST进行overwrite处理
逐一解决！首先看怎么找个合适的地方把文件送进来。这种任务当然在模式initial的时候干。
以read为关键词搜索，锁定文件IO位置。发现不少东西，比如pop设置了专门的forcing模块用于外强迫驱动。但因为forcing模块只有在单独海洋的时候才调用,为保证代码具有较高的可复用性，安全起见，我们还是首先考虑在initial module找读取外强迫场的语句。 然后呢，在initial.F90发现了这么一段可爱的内容：
!-----------------------------------------------------------------------
!
! read full 3-d t,s from input file
!
!-----------------------------------------------------------------------
case ('ccsm_startup', 'file')
first_step = .true.
if (my_task == master_task) then
write(stdout,'(a31,a)') 'Initial 3-d T,S read from file:', &
trim(init_ts_file)
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
allocate(TEMP_DATA(nx_block,ny_block,km,max_blocks_clinic))
in_file = construct_file(init_ts_file_fmt, &
full_name=trim(init_ts_file), &
record_length = rec_type_dbl, &
recl_words=nx_global*ny_global)
call data_set(in_file,'open_read')
i_dim = construct_io_dim('i',nx_global)
j_dim = construct_io_dim('j',ny_global)
k_dim = construct_io_dim('k',km)
io_temp = construct_io_field('TEMPERATURE', &
dim1=i_dim, dim2=j_dim, dim3=k_dim, &
field_loc = field_loc_center, &
field_type = field_type_scalar, &
d3d_array=TEMP_DATA)
io_salt = construct_io_field('SALINITY', &
dim1=i_dim, dim2=j_dim, dim3=k_dim, &
field_loc = field_loc_center, &
field_type = field_type_scalar, &
d3d_array=TEMP_DATA)
call data_set(in_file,'define',io_temp)
call data_set(in_file,'define',io_salt)
call data_set(in_file,'read' ,io_temp)
do iblock=1,nblocks_clinic
TRACER(:,:,:,1,curtime,iblock) = TEMP_DATA(:,:,:,iblock)
end do
call data_set(in_file,'read' ,io_salt)
do iblock=1,nblocks_clinic
TRACER(:,:,:,2,curtime,iblock) = TEMP_DATA(:,:,:,iblock)
end do
call destroy_io_field(io_temp)
call destroy_io_field(io_salt)
deallocate(TEMP_DATA)
call data_set(in_file,'close')
call destroy_file(in_file)
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,'(a12,a)') ' file read: ', trim(init_ts_file)
call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
endif
嘿嘿，和HC师兄的程序对比，如出一辙。那么我们仅仅需要在某一个合适的地方给出读取forcing场的语句就可以了。注意initial读取的时候是没有时间帧的，所以我们需要读入时间帧，并且把变量传参到step_mod或者main driver程序中。 接下来我们看如何解决这两个问题。时间问题容易解决，HC师兄的程序里给出了，直接通过dim来指定时间维即可，关键是construct_io_field这个语句，查下这个语句的定义。在io_types.F90:
function construct_io_field ( &
short_name, &
dim1, dim2, &
dim3, &
time_dim, &
long_name, &
units, &
coordinates, &
grid_loc, &
valid_range, &
field_loc, &
field_id, &
field_type, &
i0d_array, &
i1d_array, &
i2d_array, &
i3d_array, &
r0d_array, &
r1d_array, &
r2d_array, &
r3d_array, &
d0d_array, &
d1d_array, &
d2d_array, &
d3d_array) &
result (descriptor)
dim123时间都给考虑好了，意不意外，惊不惊喜？依葫芦画瓢往里送就行了。
几个神奇的参数：
character(4), intent(in), optional :: &
grid_loc ! position of field in staggered grid
integer (i4), intent(in), optional :: & ! for ghost cell updates
field_loc, &! staggering location
field_type, &! field type (scalar,vector,angle)
field_id ! previously defined id
首先的神奇之处是grid_loc，这什么鬼，虽然模板里给了3111代号，但是毕竟好奇，搜索下发现tavg.F90有如下语句：
if (present(grid_loc)) then
!*** try to decode field location from grid_loc
if (grid_loc(2:2) == '1' .and. grid_loc(3:3) == '1') then
tavg_field%field_loc = field_loc_center
else if (grid_loc(2:2) == '2' .and. grid_loc(3:3) == '2') then
tavg_field%field_loc = field_loc_NEcorner
else if (grid_loc(2:2) == '1' .and. grid_loc(3:3) == '2') then
tavg_field%field_loc = field_loc_Nface
else if (grid_loc(2:2) == '2' .and. grid_loc(3:3) == '1') then
tavg_field%field_loc = field_loc_Eface
endif
endif
东西南北？看起来是位于格点什么位置的指示。还有那些r2d_array, d2d_array代表的是数据精度real,double之类（居然没写注释！也是醉了） 查pop log发现了这个东西
2 30
________ 31
y ^ | | 32
| 3| ijk |1 33
+---> | | 34
x |________| 35
4 36
居然有这么可爱的图示……
Updated 2018-01-29
下面我们测试在main driver函数ocn_comp_mct.F90里读入一个既存的NC文件。为了效率，当然不能在run的时候每一步都读入，仅仅在initial的subroutine读入。initial和run都在ocn_mct的大module中，因此变量彼此应当可见互用。
写的时候发现一堆全局变量满天飞，POP的软件工程做得真是醉醉的……
增加namelist变量需要让bld namelist程序认识，根据这个帖子操作
增加两个subroutine，一个负责在initial的时候读取forcing文件，另一个负责在pop积分步后给出reset。
写好后发现一个大坑，由于IO module默认没有对3D var+时间的变量进行处理，如果手动增加这一功能需要改动非常多的模块，只能退而求其次，仅仅对第一层temp进行处理了。具体代码：
! ***LZN: MOD START***
subroutine ptempforcing_read
integer (i4) :: nml_error
integer (i4) :: decoupling_days=365
character (char_len) :: &
ptempf_file_name, &
ptempf_file_fmt
type (datafile) :: &
ptempf_file
type (io_field_desc) :: &
io_temp, io_ce
type (io_dim) :: &
t_dim, i_dim, j_dim, k_dim
namelist /ptempf_file_nml/ ptempf_file_name, ptempf_file_fmt
allocate(TEMP_DATA(nx_block, ny_block, decoupling_days, max_blocks_clinic))
allocate(WGT_DATA(nx_block, ny_block, decoupling_days, max_blocks_clinic))
WGT_DATA=1.0
! Read forcing file path from namelist file
if (my_task == master_task) then
open (nml_in, file=nml_filename, status='old', iostat=nml_error)
if (nml_error /= 0) then
nml_error = -1
else
nml_error = 1
endif
do while (nml_error > 0)
read(nml_in, nml=ptempf_file_nml,iostat=nml_error)
end do
if (nml_error == 0) close(nml_in)
end if
! broadcast the scalar to all MPI processors
call broadcast_scalar(ptempf_file_name, master_task)
call broadcast_scalar(ptempf_file_fmt, master_task)
ptempf_file = construct_file(ptempf_file_fmt, &
full_name=trim(ptempf_file_name),&
record_length=rec_type_real,&
recl_words=nx_global*ny_global)
call data_set(ptempf_file, 'open_read')
i_dim = construct_io_dim('i',nx_global)
j_dim = construct_io_dim('j',ny_global)
t_dim = construct_io_dim('time',decoupling_days)
io_temp = construct_io_field('TEMP_365',dim1=i_dim,dim2=j_dim,dim3=t_dim,&
long_name='Potential Temperature', units='degC', grid_loc='3111',&
field_loc = field_loc_center, &
field_type = field_type_scalar, r3d_array = TEMP_DATA)
call data_set(ptempf_file, 'define', io_temp)
call data_set(ptempf_file, 'read', io_temp)
! ru_data = construct_io_field('ru', dim1=i_dim, dim2=j_dim,&
! long_name='nudging coefficients', units='1',&
! grid_loc=' 3111', field_loc = field_loc_center,&
! field_type = field_type_scalar, r2d_array = ru)
! call data_set(sstf_file, 'define', ru_data)
! call data_set(sstf_file, 'read', ru_data)
call data_set(ptempf_file, 'close')
call destroy_io_field(io_temp)
! call destroy_io_field(ru_data)
call destroy_file(ptempf_file)
end subroutine ptempforcing_read
! fixing the SST!
subroutine ptempforcing_fix
integer (int_kind) :: iblock
! debug info
!write(stdout,*) iday_of_year, TEMP_DATA(50,50,iday_of_year,1)
! Operation: override the runtime tracer(temp)
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock = 1, nblocks_clinic
TRACER(:,:,1,1,curtime,iblock) = &
(TEMP_DATA(:,:,iday_of_year,iblock)*WGT_DATA(:,:,iday_of_year, iblock)- &
TRACER(:,:,1,1,curtime,iblock)*(1-WGT_DATA(:,:,iday_of_year, iblock)))
end do
!$OMP END PARALLEL DO
end subroutine ptempforcing_fix
! ***LZN: MOD END***
经测试一个月积分没有问题。
Updated 2018-01-29 Updated 2020-07-31