Correct allocation and assignment error for 'hgt_matrix_half' field in cospsimulator_intr.F90

[jshaw35]

Closes #1354
@cacraigucar

[brian-eaton]

This is very confusing. It appears that cosp_change_vertical_grid does want the half-level heights at the bottom of each layer. But there is a lot of code still that assumes hgt_matrix_half is dimensioned (npoints,nlevels+1). For example in function COSP_SIMULATOR we find

       rttovIN%h_surf     => cospgridIN%hgt_matrix_half(:,cospIN%Nlevels+1)

And in subroutine lidar_column we find

    real(wp),intent(in),dimension(npoints,nlevels+1) :: &
         zlev_half   ! Model half levels

I think having hgt_matrix_half only contain the bottom half levels is confusing because you don't know without searching code whether the array contains the top or the bottom levels. It's also confusing because it's unexpected. I expect the heights in hgt_matrix_half to correspond to the pressures in phalf.

My suggestion is that hgt_matrix_half be left as it is, and that the calls to cosp_change_vertical_grid be modified to pass the correct section containing the bottom half levels, e.g.,

          call cosp_change_vertical_grid(cloudsatIN%Npoints,1,cloudsatIN%Nlevels,        &
               cospgridIN%hgt_matrix(:,cloudsatIN%Nlevels:1:-1),                         &
               cospgridIN%hgt_matrix_half(:,cloudsatIN%Nlevels+1:2:-1),betamol_in,       &

[jshaw35]

@brian-eaton thanks for digging into this! I had been comparing with the current COSP version (https://github.com/CFMIP/COSPv2.0) rather than the tag that is currently used in CESM2. hgt_matrix_half is consistently dimensioned (npoints,nlevels+1) there and the calls to lidar_column and cosp_change_vertical_grid do not use the height at the top of the model.

So it seems that the changes I requested require a more recent version of COSP. @cacraigucar does it then make the most sense to wait to merge this PR until we have an official COSP tag with pending updates that do not involve RTTOV? Then this PR can include a change in the git-fleximod settings to point to an appropriate COSP tag that works with these changes. So next steps would be:

1. Merge pending changes to COSP that do not involve RTTOV. Produce a tag for CESM.
2. Point to that tag in this PR, merge.
3. Produce regression test values.
4. Test COSP-RTTOV PR and evaluate against regression tests.

[cacraigucar]

Yes, that sounds like a plan. We will need to proceed cautiously before we move forward with step 2-3 where we commit the new COSP (with the bug fix). It probably will need to have some testing by @cecilehannay before we bring it in. COSP is used at times for the CESM3 science runs, and we need to make sure it is working properly before we commit it.

[brian-eaton]

Hi @jshaw35. I've taken a look at the current head of the master branch of COSPv2.0 and see that the dimension of hgt_matrix_half has been changed to (npoints,nlevels) and that it now contains just the bottom interfaces. I agree that the path forward is to bring in an updated COSP tag which will require the changes in this PR along with appropriate changes to how hgt_matrix_half is initialized from the local variable zint.

Note also that cospsimulator_intr.F90 has an option to output the variables that are COSP inputs to the cosp_histfile_aux stream. The addfld call for ZLEV_HALF_COSP needs the vertical dimension updated.

[jshaw35]

Note for self. Relevant changes between the CESM tag and the master branch of COSP are here:
CFMIP/COSPv2.0@master...CESM_v2.1.4
CFMIP/COSPv2.0@CESM_v2.1.4...CFMIP:COSPv2.0:master
I will need to update the cospsimulator_intr.F90 to allocate and assign both hgt_matrix_half and phalf without the model top level when we use a more current tag of COSP.

[brianpm]

I agree that bringing in the updated version of COSP is the right move. I think it will also allow us to add some additional diagnostics (which I will try to look into separately).

Thanks @jshaw35 and @brian-eaton for moving this forward.

[cacraigucar]

@jshaw35 - Just a reminder that as soon as you add the COSP external update (and you have tested it with your mods), we can have the reviewers test and approve this PR and we can move it forward.

[jshaw35]

@cacraigucar I have just updated .gitmodules to reference the new COSP tag (most current COSP with no RTTOV). We will need to update the COSP driver following the discussion above.

[cacraigucar]

@cacraigucar I have just updated .gitmodules to reference the new COSP tag (most current COSP with no RTTOV). We will need to update the COSP driver following the discussion above.

@jshaw35 I assume that you will be making the change required?

[jshaw35]

@cacraigucar I have just updated .gitmodules to reference the new COSP tag (most current COSP with no RTTOV). We will need to update the COSP driver following the discussion above.

@jshaw35 I assume that you will be making the change required?

@cacraigucar yes! I don't expect it to take too long but I am sure yet when I will have the time. Hopefully sooner than later but almost certainly by September 15.

@brian-eaton the new tag for COSP should now align with the changes in this PR, resolving your points. I will review any additional changes needed to the driver since the cosp2.1.4 tag and then request an additional review.

[brianpm]

@jshaw35 -- If you have a chance, could you also see if the MODIS cloud phase diagnostics are included in the updated version of COSP? I think they should be, but we might need to add them to the CAM outputs. See, e.g., Duran et al. 2025.

[jshaw35]

@jshaw35 -- If you have a chance, could you also see if the MODIS cloud phase diagnostics are included in the updated version of COSP? I think they should be, but we might need to add them to the CAM outputs. See, e.g., Duran et al. 2025.

@brianpm , yes I don't think this should be a problem. I've just put in the new warm rain diagnostics and can add the joint histograms you're requesting. To be clear, you want the "modis_LWP_vs_ReffLIQ" and "modis_IWP_vs_ReffICE" variables, anything else?

[jshaw35]

@brianpm I've just added the MODIS LWP-REFF histograms in addition to the warm rain diagnostics from Michibata (2019). The fields are not modified by the CESM coupler (besides separating variables, naming, etc) and I've tried to follow existing naming conventions and supply units/bounds appropriately. Would you mind reviewing the output?

Average fields from a month of CAM can be found here: /glade/work/jonahshaw/COSP_RTTOV_files/COSP2.1.8_tests/cosp2.1.8_newvars.nc

[brianpm]

Thanks @jshaw35 !

Looking at the MODIS simulator, it looks like they have histograms of:

• tau-CTP ✅ (already in CAM)
• tau-CTP for liquid clouds ❓(not sure of status)
• tau-CTP for ice clouds ❓
• LWP-ReffLiq ☑︎ (you have added)
• IWP-ReffIce ☑︎
• tau-ReffIce ❓
• tau-ReffLiq ❓

I think we should add all of them if we can. What do you think?

[jshaw35]

Thanks @jshaw35 !

Looking at the MODIS simulator, it looks like they have histograms of:

* tau-CTP ✅ (already in CAM)

* tau-CTP for liquid clouds ❓(not sure of status)

* tau-CTP for ice clouds ❓

* LWP-ReffLiq ☑︎ (you have added)

* IWP-ReffIce ☑︎

* tau-ReffIce ❓

* tau-ReffLiq ❓

I think we should add all of them if we can. What do you think?

@brianpm easy to do it now! I believe those additional fields are internally named in COSP as:
modis_Optical_Thickness_vs_ReffICE (already CLRIMODIS in CESM)
modis_Optical_Thickness_vs_ReffLIQ (already CLRLMODIS in CESM)
modis_Optical_Thickness_vs_Cloud_Top_Pressure_Liq
modis_Optical_Thickness_vs_Cloud_Top_Pressure_Ice

If this looks right to you I will go ahead and add the last two histograms. Any preferences on how to name the last two fields?

[brianpm]

Yeah, that's exactly what I was looking at! Thanks!!

[jshaw35]

@brianpm all MODIS cloud histograms should now be available as output

• tau-CTP (CLMODIS)

• tau-CTP for liquid clouds (CLMODIS_LIQ)

• tau-CTP for ice clouds (CLMODIS_ICE)

• LWP-ReffLiq (IWP_REFFCLI_MODIS)

• IWP-ReffIce (LWP_REFFCLW_MODIS)

• tau-ReffIce (CLRIMODIS)

• tau-ReffLiq (CLRLMODIS)

Would you mind reviewing the outputs? Average fields from a month of CAM can be found here: /glade/work/jonahshaw/COSP_RTTOV_files/COSP2.1.8_tests/cosp2.1.8_newvars2.nc

[brianpm]

@jshaw35 -- Looking through your file, two small things.

First I think the units for the coordinate variables might be wrong:

• cosp_prs units = hPa, but values are in Pa
• cosp_refliq and cosp_refice units say microns, but it looks like the values are probably meters (?)
• cosp_iwp_modis units = "kg m-2" but values go up to 10.5, so I'm guessing this is "g m-2"?
• cosp_lwp_modis units = "kg m-2" ... also goes up to 10.125... maybe that's correct, but seems really high so maybe it is also "g m-2"

Second, I don't see CLRIMODIS and CLRLMODIS in the file.

[jshaw35]

@brianpm thanks for reviewing! Comments in-line below.

First I think the units for the coordinate variables might be wrong:

* `cosp_prs` units = hPa, but values are in Pa

Yep, you're totally right here. Now that I think about it, this has been in the CESM code for years.

* `cosp_refliq` and `cosp_refice` units say microns, but it looks like the values are probably meters (?)

Yep again!

* `cosp_iwp_modis` units = "kg m-2" but values go up to 10.5, so I'm guessing this is "g m-2"?

* `cosp_lwp_modis` units = "kg m-2" ... also goes up to 10.125... maybe that's correct, but seems really high so maybe it is also "g m-2"

These are actually correct. The upper bound on the LWP and IWP is 20 kg m-2 for both liquid and ice, so the bin center appears to be abnormally high. I think this is just meant to include all large values. The bounds arrays from cosp_config.F90 make things a bit clearer:

    real(wp),parameter,dimension(nLWP+1) :: &
         LWP_binBounds = (/0., 0.01, 0.03, 0.06, 0.10, 0.15, 0.25, 20.0/) ! kg/m2
    real(wp),parameter,dimension(nIWP+1) :: &
         IWP_binBounds = (/0., 0.02, 0.05, 0.10, 0.20, 0.40, 1.00, 20.0/) ! kg/m2

Second, I don't see CLRIMODIS and CLRLMODIS in the file.

I've updated the file with them now from the same run (so units will still be off).