Restart run output data values are different from the previous run

[ganbaranaito]

Hello everyone.

Recently, I wanted to restore SST to the climatology in the specified region. So, I modified the code 'forcing_coupled.F90' in POP source code. After modification, the model can run without any errors and the SST can be restored correctly. Unfortunately ,due to some casual issues, the run was killed in the model time 0007-02-15. Then, I switched 'CONTINUE_RUN' to 'TRUE' to restart this run from restart file 0007-01-01. However, in the restart run, I found the output data values are different from the previous run which is not killed during model time 0007-01 (even if the differences are marginal). I don't know why the output data values are different. Can anyone give me some advice? Thank you in advance!

In the code 'forcing_coupled.F90', I modified the codes after calculating STF(:,:,i,iblock):

!$OMP PARALLEL DO PRIVATE(iblock)
   do iblock = 1, nblocks_clinic
      STF(:,:,1,iblock) = (EVAP_F(:,:,iblock)*latent_heat_vapor_mks         &
                           + SENH_F(:,:,iblock) + LWUP_F(:,:,iblock)        &
                           + LWDN_F(:,:,iblock) + MELTH_F(:,:,iblock)       &
                           -(SNOW_F(:,:,iblock)+IOFF_F(:,:,iblock)) * latent_heat_fusion_mks)*  &
                             RCALCT(:,:,iblock)*hflux_factor
   enddo
 
! The following codes are which I modified:

   real (r8), dimension(1:320,1:384,1:365) ::  &
      SST_RES
   real (r8), dimension(1:320,1:384) ::  &
      lat_RES
   real (r8), dimension(1:320,1:384) ::  &
      lon_RES
   real (r8) ::  &
      tlatd,     &
      tlond,     &
      relfac,    &
      relfac1,    &
      relfac2,    &     
      tau_RES,    &
      Tm_POP
   real (r8) :: &
      Tr_RES
   integer (int_kind) ::  &
      ilon,        &
      jlat,        &
      status_RES,  &
      ncid_RES,    &
      varid_RES,   &
      i_RES,       &
      j_RES

   tau_RES = 10.*24.*3600.
   status_RES = nf90_open("/scratch/gaba/CESM/CESM2.0/CESM2_Out/clm100yrsmni.nc",nf90_nowrite,ncid_RES)
   status_RES = nf90_inq_varid(ncid_RES,"sst",varid_RES)
   status_RES = nf90_get_var(ncid_RES,varid_RES,SST_RES)
   status_RES = nf90_inq_varid(ncid_RES,"lat",varid_RES)
   status_RES = nf90_get_var(ncid_RES,varid_RES,lat_RES)
   status_RES = nf90_inq_varid(ncid_RES,"lon",varid_RES)
   status_RES = nf90_get_var(ncid_RES,varid_RES,lon_RES)
   status_RES = nf90_close(ncid_RES)
      do iblock=1,nblocks_clinic
      do jlat=1,ny_block
      do ilon=1,nx_block
        tlatd = TLAT(ilon,jlat,iblock)*radian
        tlond = TLON(ilon,jlat,iblock)*radian
         if (tlatd.gt.-20. .and. tlatd.lt.20. .and. tlond.gt.40. .and. tlond.lt.100.) then
           relfac1 = 1.0
           relfac2 = 1.0
           if (tlatd.lt.-15. .or. tlatd.gt.15.) then
             relfac1 = ABS((ABS(tlatd)-20.))*0.2
           endif
           if (tlond.lt.45. .or. tlond.gt.95.) then
             relfac2 = ABS((ABS(tlond-70.)-30.))*0.2
           endif
           relfac = MIN(relfac1,relfac2)
           do i_RES=1,320  ! lon        
                 do j_RES=1,384  ! lat            
                    if ( lat_RES(i_RES,j_RES).eq.tlatd .and. lon_RES(i_RES,j_RES).eq.tlond ) then              
                       Tr_RES = SST_RES(i_RES,j_RES,iday_of_year)
                       Tm_POP = TRACER(ilon,jlat,1,1,curtime,iblock)
                       STF(ilon,jlat,1,iblock) = STF(ilon,jlat,1,iblock) + relfac*(Tr_RES+1.-Tm_POP)*HMXL(ilon,jlat,iblock)/tau_RES                      
                    endif
                 end do
           end do
         endif
      end do
      end do
      end do

   !$OMP END PARALLEL DO

 

[klindsay]

I hypothesize that the restart failure is because you are using HMXL before it is computed in the time loop. So you are using a lagged version of HMXL. HMXL is not written to the restart file, so you're getting different HMXL on restart. You can test this hypothesis by replacing HMXL in your mods with a constant, say 50.0e2 (50 m in cm), and checking if the restart failure is fixed. Please let me know if that fixes exact restarts. If it does, then we can discuss how to include HMXL dependence and keep exact restarts working.

 

[ganbaranaito]

I hypothesize that the restart failure is because you are using HMXL before it is computed in the time loop. So you are using a lagged version of HMXL. HMXL is not written to the restart file, so you're getting different HMXL on restart. You can test this hypothesis by replacing HMXL in your mods with a constant, say 50.0e2 (50 m in cm), and checking if the restart failure is fixed. Please let me know if that fixes exact restarts. If it does, then we can discuss how to include HMXL dependence and keep exact restarts working.

Thanks for very helpful suggestions! I have replaced HMXL with a constant and found the failure is fixed! It actually disturbed me for long time, thanks! I think the next thing I should do is to include HMXL in the restart file, can you give me some tips? Thank you in advance!

Besides, in some articles' restoring experiments, HMXL tends to be fixed at 50m because 50m is the typical depth of the ocean mixed layer. So, this method may also satisfy my restoring requirement. Can I fix HMXL at 50m to do restoring experiments?

 

[klindsay]

I'm not familiar enough with the science of your restoring experiments to know if using a constant 50m HMXL is appropriate.

That said, I've implemented adding HMXL to the restart file. Modified source code, based on cesm2.1.3, is in the attached tar file. (In case you are not familiar with tar files, run 'tar xf add_HMXL_to_rest.tar' to extract the tar file contents.) In case you are not using cesm2.1.3, the tar file also has diffs of the modified files against cesm2.1.3, enabling you to apply the diffs to the version of the model that you are using.

I have verified that the attached mods build and run, add HMXL to the restart file, and read HMXL from the restart file if it is present. If HMXL is not present in the restart file, the models proceeds with HMXL=0, up to the time during the model's execution where HMXL is first computed. The enables you to run a branch or hybrid run off using a RUN_REFCASE that did not write HMXL to its restart file. If you plan to use these mods and include the HMXL dependence in your restoring formulation, please confirm that the model restarts exactly using the mods.

 

[ganbaranaito]

I'm not familiar enough with the science of your restoring experiments to know if using a constant 50m HMXL is appropriate.

That said, I've implemented adding HMXL to the restart file. Modified source code, based on cesm2.1.3, is in the attached tar file. (In case you are not familiar with tar files, run 'tar xf add_HMXL_to_rest.tar' to extract the tar file contents.) In case you are not using cesm2.1.3, the tar file also has diffs of the modified files against cesm2.1.3, enabling you to apply the diffs to the version of the model that you are using.

I have verified that the attached mods build and run, add HMXL to the restart file, and read HMXL from the restart file if it is present. If HMXL is not present in the restart file, the models proceeds with HMXL=0, up to the time during the model's execution where HMXL is first computed. The enables you to run a branch or hybrid run off using a RUN_REFCASE that did not write HMXL to its restart file. If you plan to use these mods and include the HMXL dependence in your restoring formulation, please confirm that the model restarts exactly using the mods.

Thank you so much! The version of CESM which I use is also CESM2.1.3. I will give it a try as soon as possible. Thanks again!

 

[ganbaranaito]

I'm not familiar enough with the science of your restoring experiments to know if using a constant 50m HMXL is appropriate.

That said, I've implemented adding HMXL to the restart file. Modified source code, based on cesm2.1.3, is in the attached tar file. (In case you are not familiar with tar files, run 'tar xf add_HMXL_to_rest.tar' to extract the tar file contents.) In case you are not using cesm2.1.3, the tar file also has diffs of the modified files against cesm2.1.3, enabling you to apply the diffs to the version of the model that you are using.

I have verified that the attached mods build and run, add HMXL to the restart file, and read HMXL from the restart file if it is present. If HMXL is not present in the restart file, the models proceeds with HMXL=0, up to the time during the model's execution where HMXL is first computed. The enables you to run a branch or hybrid run off using a RUN_REFCASE that did not write HMXL to its restart file. If you plan to use these mods and include the HMXL dependence in your restoring formulation, please confirm that the model restarts exactly using the mods.

Hello, Keith!
I have used your modified codes in my run and found HMXL is written into the restart files! The restart run output data values are same with the previous run! Your modified codes solved my problems completely. Thanks a lot!

 

[klindsay]

Thanks for the update. I'm glad to hear that it works for you.

In case someone reads this thread later and attempts to add another variable to the restart file, I'm including some details of the implementation of writing and reading HMXL to and from the restart file. There are some steps in doing this where it is easy to make a mistake or get stuck.
1) HMXL is only allocated and computed by the ocean model if vmix_choice=='kpp'. So the mods only write and read HMXL to and from the restart file if that is the case. The mods check this by checking if HMXL is allocated.
2) Adding HMXL to write_restart in restart.F90, inside a check for allocated(HMXL) is straightforward. The mods follow the example of QFLUX, which is also conditionally written to the restart file. The module restart uses the module vmix_kpp to access its module variable HMXL.
3) Adding HMXL to read_restart in restart.F90, with a check for allocated(HMXL) doesn't work, because HMXL is allocated during initialization after read_restart is called. The allocation status of HMXL is not known when read_restart is called. So I added the read for HMXL to init_vmix_kpp in vmix_kpp.F90, after HMXL is allocated.
4) HMXL should only be read from the restart file if the model is initializing from a restart (i.e., a branch or hybrid run, not a startup run). The code mods check this by checking to see if the input restart filename has been set.
5) vmix_kpp.F90 cannot use the module restart to access the input restart filename and format module variables, because this would lead to the modules restart and vmix_kpp using each other, and this circular dependence causes problems for building the model. So the restart filename and format are passed down the subroutine call tree from pop_init_phase2 in initial.F90 to init_vertical_mix in vertical_mix.F90 to init_vertical_mix in vmix_kpp.F90. Note that the input restart file format is init_ts_file_fmt, not restart_fmt. This latter variable is used for the output restart file format, which has the potential to differ from the input restart file format.