Question about modifying ocean surface albedo in SourceCode

[anserwei]

Dear CESM team,

Hello, I want to modify the source code of ocean surface albedo using a new calculation method for one of my Ph.D. research topics.

In CESM 2.1.3, the source code of ocean surface albedo is in the subroutine seq_flux_ocnalb_mct( infodata, ocn, a2x_o, fractions_o, xao_o ) of seq_flux_mct.F90 (source code path: ./cime/src/drivers/mct/main/seq_flux_mct.F90):

There are two chooses for calculating the ocean surface albedo:
1. one is daily mean without considering the zenith angle dependence (line 821-843 in seq_flux_mct.F90):

if (flux_albav) then do n=1,nloc_o anidr = seq_flux_mct_albdir avsdr = seq_flux_mct_albdir anidf = seq_flux_mct_albdif avsdf = seq_flux_mct_albdif ! Albedo is now function of latitude (will be new implementation) !rlat = const_deg2rad * lats(n) !anidr = 0.069_r8 - 0.011_r8 * cos(2._r8 * rlat) !avsdr = anidr !anidf = anidr !avsdf = anidr xao_o%rAttr(index_xao_So_avsdr,n) = avsdr xao_o%rAttr(index_xao_So_anidr,n) = anidr xao_o%rAttr(index_xao_So_avsdf,n) = avsdf xao_o%rAttr(index_xao_So_anidf,n) = anidf end do update_alb = .true. else

2. the other one is instantaneous with the zenith angle dependence (line 876 - 901 in seq_flux_mct.F90):

! Compute albedos do n=1,nloc_o rlat = const_deg2rad * lats(n) rlon = const_deg2rad * lons(n) cosz = shr_orb_cosz( calday, rlat, rlon, delta ) if (cosz > 0.0_r8) then !--- sun hit -- anidr = (.026_r8/(cosz**1.7_r8 + 0.065_r8)) + & (.150_r8*(cosz - 0.100_r8 ) * & (cosz - 0.500_r8 ) * & (cosz - 1.000_r8 ) ) avsdr = anidr anidf = seq_flux_mct_albdif avsdf = seq_flux_mct_albdif else !--- dark side of earth --- anidr = 1.0_r8 avsdr = 1.0_r8 anidf = 1.0_r8 avsdf = 1.0_r8 end if xao_o%rAttr(index_xao_So_avsdr,n) = avsdr xao_o%rAttr(index_xao_So_anidr,n) = anidr xao_o%rAttr(index_xao_So_avsdf,n) = avsdf xao_o%rAttr(index_xao_So_anidf,n) = anidf end do ! nloc_o


By default, the Coupler uses the second one, i.e., daily mean downward solar radiation * instantaneous surface albedo. On the dark side, there is no solar absorption, so the net solar radiation send to the ocean components is much smaller than that expected.

I write a simple subroutine that ocean surface albedo is parameterized as function of four variables:

Inputs:
Cosine of solar zenith angle: cosz
ocean surface chlorophyll concentration: chl (unit: mg/m-3); total chlorophyll
0meter wind speed: wspd10 (unit: m/s)；
wind direction:wspd_dir (unit: degree);
wavelength start: wl_s (unit: nm);
wavelength end: wl_e (unit: nm);

outputs:
direct albedo: adr
diffuse albedo: adf

The subroutine looks like this: call osa_para(cosz, chl, wspd10, wspd_dir, wl_s, wl_e, adr, adf)

I want to replace the second option of ocean surface albedo calculation with my own subroutine osa_para,
(i.e., replaced

anidr = (.026_r8/(cosz**1.7_r8 + 0.065_r8)) + & (.150_r8*(cosz - 0.100_r8 ) * & (cosz - 0.500_r8 ) * & (cosz - 1.000_r8 ) ) avsdr = anidr anidf = seq_flux_mct_albdif avsdf = seq_flux_mct_albdif
)

For instance, with the subroutine osa_para,

all osa_para(cosz, chl, wspd10, wspd_dir, 400, 700, avsdr, avsdf)

Parameters/availabilities for calling the osa_para(cosz, chl, wspd10, wspd_dir, wl_s, wl_e, adr, adf):

(1). ‘Cosz’ is provided directly in .//cime/src/drivers/mct/main/seq_flux_mct.F90.
(2). ‘chl’ is not available in .//cime/src/drivers/mct/main/seq_flux_mct.F90.
(3). ‘wspd10’ is the variable ‘duu10n’ in .//cime/src/drivers/mct/main/seq_flux_mct.F90.
(4). ‘wspd_dir’ can be calculated from varables ‘ubot ’ and ‘vbot’ in .//cime/src/drivers/mct/main/seq_flux_mct.F90.
(5) ‘wl_s’ and ‘wl_s’ are fixed integer values. (e.g., wl_s = 300, wl_e = 700)


Questions:

(1). How to add the ocean surface chlorophyll concentration ‘chl’ to the seq_flux_mct.F90 file (specifically, in the subroutine seq_flux_ocnalb_mct( infodata, ocn, a2x_o, fractions_o, xao_o))?

The subroutine seq_flux_ocnalb_mct in the seq_flux_mct.F90 is written as:
!=============================================================================== subroutine seq_flux_ocnalb_mct( infodata, ocn, a2x_o, fractions_o, xao_o ) !----------------------------------------------------------------------- ! Arguments type(seq_infodata_type) , intent(in) :: infodata type(component_type) , intent(in) :: ocn type(mct_aVect) , intent(in) :: a2x_o type(mct_aVect) , intent(inout) :: fractions_o type(mct_aVect) , intent(inout) :: xao_o

......

If ‘chl’ is added to this subroutine seq_flux_ocnalb_mct. The input parameters: infodata, ocn and a2x_o, need to be changed ?? how?

(2) If the 'chl' is provided by the pop, what additional process is needed (like CN_CHL_TYPE=diagnostic??? )?
What about if the monthly mean 'chl' are from Oceancolor satellite observation climatology dataset? what additional process is needed?

(3) So, the namlist of CIME is needed to be changed accordingly? how?

It will be very appreciated if anyone could guide me on how to implement modifications (mentioned above).

Thank you very much!


Best,
Johnny Wong

 

[erik]

Since you are talking about chlorophyll concentration in coupler ocean code, I think this would be better served in the Ocean forum about MARBL I think what you'll have to do is pass that concentration from MARBL and the ocean model to the coupler code. But, ocean modelers are going to be better at helping you through this.

 

[mlevy]

Hi Johnny,

This is a very complicated process -- I'll outline it below, and provide details where I know them, but we may need to reach out to other people who are more familiar with coupler code than I am to fill in the gaps. Basically, the steps you'll need to take are

1. Tell the coupler that the ocean will be passing chlorophyll (I don't know exactly how to do this, but maybe look for 'So_u' in the coupler code to see how it handles recieving POP's UVEL variable?)
2. Tell POP to provide chlorophyll in the buffer passed to the coupler

* This occurs in the POP_sum_buffer() subroutine in drivers/cpl/ocn_import_export.F90.
* Specifically, after step (1), the third dimension of the array SBUFF_SUM will be longer and you will need to add another entry in the loop that updates this array. In the sample code below, I'm assuming you modified the coupler

   SBUFF_SUM(:,:,iblock,index_o2x_So_chl) =   &
      SBUFF_SUM(:,:,iblock,index_o2x_So_chl) + delt*  &
                                   CHL(:,:,curtime,iblock)

* Two notes: (1) you'll need to include a use sw_absorption, only : CHL statement at the top of the module (and give CHL the public attribute in that module), and (2) you'll need to define index_o2x_So_chl in drivers/cpl/POP_CplIndices.F90, with something like

index_o2x_So_chl        = mct_avect_indexra(o2x,'So_chl')

Where I assume the 'So_chl' needs to match something you defined in the coupler to let it know that chlorophyll will be coming.

3. At this point, you can go back to coupler and again follow So_u as an example to figure out how to store the chlorophyll in memory and access it in your subroutines.

Let me know if you run into issues - I suspect you will, as I haven't tested the above code snippets at all... but it should at least provide a reasonable guide for how to proceed.

 

[anserwei]

Hi Johnny,

This is a very complicated process -- I'll outline it below, and provide details where I know them, but we may need to reach out to other people who are more familiar with coupler code than I am to fill in the gaps. Basically, the steps you'll need to take are

1. Tell the coupler that the ocean will be passing chlorophyll (I don't know exactly how to do this, but maybe look for 'So_u' in the coupler code to see how it handles recieving POP's UVEL variable?)
2. Tell POP to provide chlorophyll in the buffer passed to the coupler

* This occurs in the POP_sum_buffer() subroutine in drivers/cpl/ocn_import_export.F90.
* Specifically, after step (1), the third dimension of the array SBUFF_SUM will be longer and you will need to add another entry in the loop that updates this array. In the sample code below, I'm assuming you modified the coupler

   SBUFF_SUM(:,:,iblock,index_o2x_So_chl) =   &
      SBUFF_SUM(:,:,iblock,index_o2x_So_chl) + delt*  &
                                   CHL(:,:,curtime,iblock)

* Two notes: (1) you'll need to include a use sw_absorption, only : CHL statement at the top of the module (and give CHL the public attribute in that module), and (2) you'll need to define index_o2x_So_chl in drivers/cpl/POP_CplIndices.F90, with something like

index_o2x_So_chl        = mct_avect_indexra(o2x,'So_chl')

Where I assume the 'So_chl' needs to match something you defined in the coupler to let it know that chlorophyll will be coming.

3. At this point, you can go back to coupler and again follow So_u as an example to figure out how to store the chlorophyll in memory and access it in your subroutines.

Let me know if you run into issues - I suspect you will, as I haven't tested the above code snippets at all... but it should at least provide a reasonable guide for how to proceed.

Hi, Michael,

Thank you for your suggestions.

The motivations of my study (20 years GCM run):

(1) Plan to provide two options of chlorophyll concentration in the 'seq_flux_mct.F90': (1-1). the diagnostic (1-2). The prescribe.
(2) Modify the ocean surface albedo calculation in CESM using the new algorithm.
(3) Compare ocean surface albedo results of the new algorithm to the original algorithm in CESM, and analyze the influence of chlorophyll concentration on ocean surface albedo.
Specifically, run CESM twice: one outputs the ocean surface albedo of the new algorithm; another outputs the ocean surface albedo of the original algorithm.
(4) Run CESM twice: access changes of atmospheric parameters between the original and the new algorithm.
For instance, atmospheric parameters include air temperature, precipitation, specific humidity, aerosol optical depth, downward and upward radiation, etc.
(5) For ‘compsets’, I don’t have a preference yet. But, I think the selected ‘compset’ should be a coupled air-ocean-seaice-land ‘compset’. The changes in atmospheric parameters should not be limited to the atmosphere over the ocean.

Currently, I am modifying the code following your suggestions. I will post some updates once I have.

Thank you very much.

 

[Riley]

hi, i am interested that whether "subroutine seq_flux_ocnalb_mct" also acts on slab ocn model or not?

 

[anserwei]

Hi Johnny,

This is a very complicated process -- I'll outline it below, and provide details where I know them, but we may need to reach out to other people who are more familiar with coupler code than I am to fill in the gaps. Basically, the steps you'll need to take are

1. Tell the coupler that the ocean will be passing chlorophyll (I don't know exactly how to do this, but maybe look for 'So_u' in the coupler code to see how it handles recieving POP's UVEL variable?)
2. Tell POP to provide chlorophyll in the buffer passed to the coupler

* This occurs in the POP_sum_buffer() subroutine in drivers/cpl/ocn_import_export.F90.
* Specifically, after step (1), the third dimension of the array SBUFF_SUM will be longer and you will need to add another entry in the loop that updates this array. In the sample code below, I'm assuming you modified the coupler

   SBUFF_SUM(:,:,iblock,index_o2x_So_chl) =   &
      SBUFF_SUM(:,:,iblock,index_o2x_So_chl) + delt*  &
                                   CHL(:,:,curtime,iblock)

* Two notes: (1) you'll need to include a use sw_absorption, only : CHL statement at the top of the module (and give CHL the public attribute in that module), and (2) you'll need to define index_o2x_So_chl in drivers/cpl/POP_CplIndices.F90, with something like

index_o2x_So_chl        = mct_avect_indexra(o2x,'So_chl')

Where I assume the 'So_chl' needs to match something you defined in the coupler to let it know that chlorophyll will be coming.

3. At this point, you can go back to coupler and again follow So_u as an example to figure out how to store the chlorophyll in memory and access it in your subroutines.

Let me know if you run into issues - I suspect you will, as I haven't tested the above code snippets at all... but it should at least provide a reasonable guide for how to proceed.

=================================================================================================================

HI, Michael

Hope you are doing well. Following your suggestions, I modified the codes in CESM version 2.1.3. I want to run a case to see if my modifications are valid or not. My case could successfully build and run with no error. However, the case was crashed at a certain point (it's around 20mins after the submission).

In the CaseStatus:
case.run error
ERROR: RUN FAIL: Command 'mpirun -np 672 /scratch/user/anser/b.e20.B1850.f19_g17.test.mdf/bld/cesm.exe >> cesm.log.$LID 2>&1 ' failed
See log file for details: /scratch/user/anser/b.e20.B1850.f19_g17.test.mdf/run/cesm.log.7997926.210511-153045

In the cesm.log... files:
ERROR:
component_mod:check_fields NaN found in ATM instance: 1 field Sa_z 1d global
......

Might you have some clue where to start for the debugging process? Could you help me fix this? Thank you very much!

I provided some info here that is relative to my case

%%%%%=====
1. Every step I took:

(1)

./create_newcase --case $SCRATCH/case_dir/b.e20.B1850.f19_g17.test.mdf --compset B1850 --res f19_g17 --mach terra --run-unsupported

(2)

cd b.e20.B1850.f19_g17.test.mdf

(3)

./xmlchange DIN_LOC_ROOT=$SCRATCH/input_data_dir/b.e20.B1850.f19_g17/

(4)

./xmlchange NTASKS_ATM=224,NTHRDS_ATM=2,ROOTPE_ATM=0
./xmlchange NTASKS_ICE=112,NTHRDS_ICE=2,ROOTPE_ICE=112
./xmlchange NTASKS_LND=112,NTHRDS_LND=2,ROOTPE_LND=0
./xmlchange NTASKS_CPL=224,NTHRDS_CPL=2,ROOTPE_CPL=0
./xmlchange NTASKS_ROF=112,NTHRDS_ROF=2,ROOTPE_ROF=0
./xmlchange NTASKS_OCN=112,NTHRDS_OCN=2,ROOTPE_OCN=224
./xmlchange NTASKS_GLC=56,NTHRDS_GLC=2,ROOTPE_GLC=0
./xmlchange NTASKS_WAV=56,NTHRDS_WAV=2,ROOTPE_WAV=0
 ./xmlchange NTASKS_ESP=1,NTHRDS_ESP=2,ROOTPE_ESP=0

(5)

./xmlchange RUN_TYPE=hybrid,RUN_STARTDATE=2010-03-01,STOP_OPTION=nmonths,STOP_N=6,DOUT_S=TRUE

(6)

./case.setup

(7)

(7-1) code modification

copy modified *.F90 to .//SourceMods/src.pop/
copy modified *.F90 to .//SourceMods/src.drv/

(7-2) edit the user_nl_cam and add the lines:

nhtfrq = 0, 0
mfilt = 1, 1
fincl2 = 'ASDIR','ASDIF','ALDIR','ALDIF'

(7-3) edit the user_nl_pop and add the lines:
tavg_freq_opt = 'nmonth'
tavg_freq = 1

(7-4) edit the user_nl_cice and add the lines:

histfreq = 'm'
histfreq_n = 1

(7-5) edit the user_nl_clm and add the lines:

hist_nhtfrq = 0
hist_mfilt = 1

(8)

./case.build --skip-provenance-check

(9)

./case.submit

%%%%%=====
2. I attach some files here that may be helpful. see attached files.

For instance, files include CaseStatus, README.case, config_compilers.xml, config_machines.xml, and config_batch.xml files, all log files.

%%%%%=====
3. I found this link that describes a similar CESM issue (the last reply):

https://bb.cgd.ucar.edu/cesm/threads/error%EF%BC%9Acomponent_mod-check_fields-nan-found-in-atm-instance-1-field-sa_z-1d-global.5797/#post-38962

But, its solution doesn't make sense to me.

%%%%%=====
4. I don't know if the error is relative to the compset I used (I used B1850):
You mentioned that If a compset that sets POP2%ECO, then MARBL will provide BGC tracers and diagnostics to POP (including chlorophyll).
So, could you provide some compsets you suggested for my studies?

looking forward to your reply. Thank you very much!

Johnny

 

[ganbaranaito]

=================================================================================================================

HI, Michael

Hope you are doing well. Following your suggestions, I modified the codes in CESM version 2.1.3. I want to run a case to see if my modifications are valid or not. My case could successfully build and run with no error. However, the case was crashed at a certain point (it's around 20mins after the submission).

In the CaseStatus:
case.run error
ERROR: RUN FAIL: Command 'mpirun -np 672 /scratch/user/anser/b.e20.B1850.f19_g17.test.mdf/bld/cesm.exe >> cesm.log.$LID 2>&1 ' failed
See log file for details: /scratch/user/anser/b.e20.B1850.f19_g17.test.mdf/run/cesm.log.7997926.210511-153045

In the cesm.log... files:
ERROR:
component_mod:check_fields NaN found in ATM instance: 1 field Sa_z 1d global
......

Might you have some clue where to start for the debugging process? Could you help me fix this? Thank you very much!

I provided some info here that is relative to my case

%%%%%=====
1. Every step I took:

(1)

./create_newcase --case $SCRATCH/case_dir/b.e20.B1850.f19_g17.test.mdf --compset B1850 --res f19_g17 --mach terra --run-unsupported

(2)

cd b.e20.B1850.f19_g17.test.mdf

(3)

./xmlchange DIN_LOC_ROOT=$SCRATCH/input_data_dir/b.e20.B1850.f19_g17/

(4)

./xmlchange NTASKS_ATM=224,NTHRDS_ATM=2,ROOTPE_ATM=0
./xmlchange NTASKS_ICE=112,NTHRDS_ICE=2,ROOTPE_ICE=112
./xmlchange NTASKS_LND=112,NTHRDS_LND=2,ROOTPE_LND=0
./xmlchange NTASKS_CPL=224,NTHRDS_CPL=2,ROOTPE_CPL=0
./xmlchange NTASKS_ROF=112,NTHRDS_ROF=2,ROOTPE_ROF=0
./xmlchange NTASKS_OCN=112,NTHRDS_OCN=2,ROOTPE_OCN=224
./xmlchange NTASKS_GLC=56,NTHRDS_GLC=2,ROOTPE_GLC=0
./xmlchange NTASKS_WAV=56,NTHRDS_WAV=2,ROOTPE_WAV=0
./xmlchange NTASKS_ESP=1,NTHRDS_ESP=2,ROOTPE_ESP=0

(5)

./xmlchange RUN_TYPE=hybrid,RUN_STARTDATE=2010-03-01,STOP_OPTION=nmonths,STOP_N=6,DOUT_S=TRUE

(6)

./case.setup

(7)

(7-1) code modification

copy modified *.F90 to .//SourceMods/src.pop/
copy modified *.F90 to .//SourceMods/src.drv/

(7-2) edit the user_nl_cam and add the lines:

nhtfrq = 0, 0
mfilt = 1, 1
fincl2 = 'ASDIR','ASDIF','ALDIR','ALDIF'

(7-3) edit the user_nl_pop and add the lines:
tavg_freq_opt = 'nmonth'
tavg_freq = 1

(7-4) edit the user_nl_cice and add the lines:

histfreq = 'm'
histfreq_n = 1

(7-5) edit the user_nl_clm and add the lines:

hist_nhtfrq = 0
hist_mfilt = 1

(8)

./case.build --skip-provenance-check

(9)

./case.submit

%%%%%=====
2. I attach some files here that may be helpful. see attached files.

For instance, files include CaseStatus, README.case, config_compilers.xml, config_machines.xml, and config_batch.xml files, all log files.

%%%%%=====
3. I found this link that describes a similar CESM issue (the last reply):

But, its solution doesn't make sense to me.

%%%%%=====
4. I don't know if the error is relative to the compset I used (I used B1850):
You mentioned that If a compset that sets POP2%ECO, then MARBL will provide BGC tracers and diagnostics to POP (including chlorophyll).
So, could you provide some compsets you suggested for my studies?

looking forward to your reply. Thank you very much!

Johnny

Hi Johnny,

Have you solved this problem? I also met this error. Can you give me some suggestions? Thank you!

Gaya