search

ufs-community / UFS_UTILS

Utilities for the NCEP models.
Other
21 stars 104 forks source link

unexplained constants in chgres_cube.fd/program_setup.F90, subroutine calc_soil_params() #374

Open edwardhartnett opened 3 years ago

edwardhartnett commented 3 years ago

In in chgres_cube.fd/program_setup.F90, subroutine calc_soil_params(), we have:

  SATDW(I)  = BB(I)*SATDK(I)*(SATPSI(I)/MAXSMC(I))
   REFSMC1 = MAXSMC(I)*(5.79E-9/SATDK(I)) **(1.0/(2.0*BB(I)+3.0))
   REFSMC(I) = REFSMC1 + (MAXSMC(I)-REFSMC1) / SMHIGH
   WLTSMC1 = MAXSMC(I) * (200.0/SATPSI(I))**(-1.0/BB(I))
   WLTSMC(I) = WLTSMC1 - SMLOW * WLTSMC1

What is 5.79E-9? 1.0? 2.0? 3.0? 200.0? -1.0?

These should be named constants. (At least some of them.)

Also some comments would be nice.

edwardhartnett commented 3 years ago

I looked in this paper to try and figure out the formulas: https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2002JD003296.

Seems like section 8.1 of the paper is the relevant part.

FX = (theta1 - theta_dry)/(theta_sat - theta_dry) and

Edir = (1 - sigma_f)(FX^^fx)(Ep)

fx = 2

I'm having a hard time matching the equations to the code. ;-)

edwardhartnett commented 3 years ago

I'm also looking at the equations on page 573 of https://journals.ametsoc.org/view/journals/mwre/129/4/1520-0493_2001_129_0569_caalsh_2.0.co_2.xml. I see some terms there also used in the code, like

hydraulic conductivity = K = satdk soil water diffusivity = D = satdw

GeorgeGayno-NOAA commented 3 years ago

@barlage Can you or someone on the land team help?

edwardhartnett commented 3 years ago

I also see these numbers:

data bb_statsgo /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
            6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
            5.25, -9.99, 4.05, 4.26/

 data maxsmc_statsgo /0.395, 0.421, 0.434, 0.476, 0.476, 0.439, &
              0.404, 0.464, 0.465, 0.406, 0.468, 0.457, &
              0.464, -9.99, 0.200, 0.421/

 data satdk_statsgo /1.7600e-4, 1.4078e-5, 5.2304e-6, 2.8089e-6, 2.8089e-6, &
             3.3770e-6, 4.4518e-6, 2.0348e-6, 2.4464e-6, 7.2199e-6, &
             1.3444e-6, 9.7384e-7, 3.3770e-6,     -9.99, 1.4078e-5, &
             1.4078e-5/

 data satpsi_statsgo /0.0350, 0.0363, 0.1413, 0.7586, 0.7586, 0.3548, &
              0.1349, 0.6166, 0.2630, 0.0977, 0.3236, 0.4677, &
              0.3548, -9.99,  0.0350, 0.0363/

 data bb_zobler /4.26,  8.72, 11.55,  4.74, 10.73,  8.17, &
                 6.77,  5.25,  4.26/

 data maxsmc_zobler /0.421, 0.464, 0.468, 0.434, 0.406, 0.465, &
                     0.404, 0.439, 0.421/

 data satdk_zobler /1.41e-5, 0.20e-5, 0.10e-5, 0.52e-5, 0.72e-5, &
                    0.25e-5, 0.45e-5, 0.34e-5, 1.41e-5/

 data satpsi_zobler /0.040, 0.620, 0.470, 0.140, 0.100, 0.260,  &
                     0.140, 0.360, 0.040/

I am wondering where they come from. I see in this paper some similar numbers, but they are not the same. Where our numbers derived from the paper's? If so, how?

https://pdfs.semanticscholar.org/e3a3/d0e56552b17c719c5a50dabc10ae471bf821.pdf

barlage commented 3 years ago

I can help out with some of this, but not immediately. I will try to respond before tomorrow.

edwardhartnett commented 3 years ago

It's not urgent. I'm writing tests for this code and just want to check the numbers, rather than accepting them as correct without question.

edwardhartnett commented 3 years ago

@barlage any progress here?

edwardhartnett commented 3 years ago

Any progress here?

edwardhartnett commented 3 years ago

@barlage any progress here? I am somewhat concerned that these numbers are so mysterious. Do we understand where they came from?

edwardhartnett commented 3 years ago

@arunchawla-NOAA perhaps you can jump in and try and shake loose some answers from the land team? There are a bunch of constants in the code, and we're trying to learn where they came from. They were attributed to some science papers, but do not appear in those papers. Some additional work was done to derive these numbers - what was it?

Perhaps we simply don't know where these numbers came from? What do we do then?

barlage commented 3 years ago

This issue seems minor relative to other issues I'm dealing with so has been superseded by other concerns. If this is holding up some larger development or truly requires management intervention, let me know.

HelinWei-NOAA commented 3 years ago

This portion of codes is copied from set_soilveg.f, which is part of Noah LSM in the UFS. I don't think you need to know the detail here. You can just refer to this paper: Ek, M. B., K. E. Mitchell, Y. Lin, E. Rogers, P. Grunmann, V. Koren, G. Gayno, and J. D. Tarpley, Implementation of Noah land surface model advances in the National Centers for Environmental Prediction operational mesoscale Eta model, J. Geophys. Res., 108(D22), 8851, doi:10.1029/2002JD003296, 2003

Make sure whenever we have updates in set_soilveg.f, you need to keep the same here.

edwardhartnett commented 3 years ago

The paper you cite has been added to the documentation. I have read the relevant section and it does not provide the numbers used in the code.

The goal here is testing, which is to make sure these numbers are correct.

Where is the file set_soilveg.f which you refer to? That is, what repository and directory?

HelinWei-NOAA commented 3 years ago

https://github.com/HelinWei-NOAA/ccpp-physics/blob/master/physics/set_soilveg.f

arunchawla-NOAA commented 3 years ago

Copying files from one place to another is not a good strategy. If there are routines for land set up that are needed in both the physics module and in chgres then we need to have a smarter strategy. Perhaps a library that can be called in both places. This does not have to be solved right now.

At the very least these constants need to be documented somewhere. Are they documented as part of set_soilveg?

HelinWei-NOAA commented 3 years ago

Those parameters are documented as part of Noah LSM in https://github.com/HelinWei-NOAA/ccpp-physics/blob/master/physics/sflx.f

edwardhartnett commented 3 years ago

OK, thanks, that is helpful.

We will ponder how to best handle constants taken from other codes.

edwardhartnett commented 3 years ago

OK, let's confine the discussion to one set of parameters, so we have a specific example.

In the chgres_cube file program_setup.f we have:

! using stasgo table
 data bb_statsgo /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
            6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
            5.25, -9.99, 4.05, 4.26/

So are those the correct numbers? @HelinWei-NOAA you provided this link to a file on your fork: https://github.com/HelinWei-NOAA/ccpp-physics/blob/master/physics/set_soilveg.f

When I look in there, I find this:

      bb         =(/4.26,  8.72, 11.55, 4.74, 10.73,  8.17,
     &            6.77,  5.25,  4.26, 0.00,  0.00,  0.00,
     &            0.00,  0.00,  0.00, 0.00,  0.00,  0.00,
     &            0.00,  0.00,  0.00, 0.00,  0.00,  0.00,
     &            0.00,  0.00,  0.00, 0.00,  0.00,  0.00/)

When I look at this file on the master branch I see:

BB         =(/4.05,  4.26, 4.74, 5.33, 5.33,  5.25,
     &            6.77,  8.72,  8.17, 10.73, 10.39,  11.55,
     &            5.25,  4.26,  4.05, 4.26,  11.55,  4.05,
     &            4.05,  0.00,  0.00, 0.00,  0.00,  0.00,
     &            0.00,  0.00,  0.00, 0.00,  0.00,  0.00/)

There is no develop branch.

None of these numbers match the others. ;-)

Some questions:

  1. Should we be looking at your fork for the authoritative numbers?
  2. Am I looking at the correct arrays?
  3. The chgres_cube code has 16 elements in this array, the set_soilveg.f90 has 30. Is that expected?
  4. What numbers should be used?
HelinWei-NOAA commented 3 years ago

Hi Edward,

There are two different soil type data which can be used in the model. The previous (old) one is zolber with 9 classifications and the other one is statsgo with 16 classifications. So in that subroutine the parameters for both types are defined.

For example,

if(isot.eq.0) then (this one is for old soil type data)

  bb         =(/4.26,  8.72, 11.55, 4.74, 10.73,  8.17,
 &            6.77,  5.25,  4.26, 0.00,  0.00,  0.00,
 &            0.00,  0.00,  0.00, 0.00,  0.00,  0.00,
 &            0.00,  0.00,  0.00, 0.00,  0.00,  0.00,
 &            0.00,  0.00,  0.00, 0.00,  0.00,  0.00/)

..... else

! using stasgo table BB =(/4.05, 4.26, 4.74, 5.33, 5.33, 5.25, & 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, & 5.25, 4.26, 4.05, 4.26, 11.55, 4.05, & 4.05, 0.00, 0.00, 0.00, 0.00, 0.00, & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/)

Here in charge_cube, we assume the statsgo data is used.

edwardhartnett commented 3 years ago

OK, thanks.

Should we be looking at a fork for these numbers?

HelinWei-NOAA commented 3 years ago

At this time, Yes. But in the future we will move this part from the code to a table file, then we should find a way to let charge_cube and UFS to share this file.

edwardhartnett commented 3 years ago

OK. In that fork I see this:

! using stasgo table
      BB         =(/4.05,  4.26, 4.74, 5.33, 5.33,  5.25,
     &            6.77,  8.72,  8.17, 10.73, 10.39,  11.55,
     &            5.25,  4.26,  4.05, 4.26,  11.55,  4.05,
     &            4.05,  0.00,  0.00, 0.00,  0.00,  0.00,
     &            0.00,  0.00,  0.00, 0.00,  0.00,  0.00/)

These are lines 245-250. Are these the correct numbers?

In chgres_cube we have:

data bb_statsgo /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
            6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
            5.25, -9.99, 4.05, 4.26/

The last four numbers are different. Is this expected?

HelinWei-NOAA commented 3 years ago

Type 14 is water, so -9,99 is used in chgres_cube. Actually statsgo soil type has 19 classifications, here we only use the first 16. George, did you do this for some reason?

! upgraded to statsgo (19-type) ! 1: sand ! 2: loamy sand ! 3: sandy loam ! 4: silt loam ! 5: silt ! 6:loam ! 7:sandy clay loam ! 8:silty clay loam ! 9:clay loam ! 10:sandy clay ! 11: silty clay ! 12: clay ! 13: organic material ! 14: water ! 15: bedrock ! 16: other (land-ice) ! 17: playa ! 18: lava ! 19: white sand

On Tue, Apr 27, 2021 at 3:21 PM Edward Hartnett @.***> wrote:

OK. In that fork I see this:

! using stasgo table BB =(/4.05, 4.26, 4.74, 5.33, 5.33, 5.25, & 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, & 5.25, 4.26, 4.05, 4.26, 11.55, 4.05, & 4.05, 0.00, 0.00, 0.00, 0.00, 0.00, & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/)

These are lines 245-250. Are these the correct numbers?

In chgres_cube we have:

data bb_statsgo /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, & 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, & 5.25, -9.99, 4.05, 4.26/

The last four numbers are different. Is this expected?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/UFS_UTILS/issues/374#issuecomment-827856365, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALPHKYHVNYL6JGGDPT4YTSDTK4FD5ANCNFSM4YPFFIWA .

GeorgeGayno-NOAA commented 3 years ago

@HelinWei-NOAA None of our input datasets use classifications 17 thru 19. And none of our standard grids in ./fix use those classifications either. So, I did not include them.

Does anyone use them? Should I add them?

HelinWei-NOAA commented 3 years ago

I am not aware that anyone uses them. But for consistence, maybe we should add them.