CAM5 MG1.0 liquid and ice number concentrations much smaller than in evaluation paper

[andreasb]

I am comparing liquid and ice number concentrations from CESM1/CAM5 runs (v1_2_1, MG1.0, GEOS5-nudged and free-running setups, present day, e.g. user_compset 2000_CAM5_CLM40%CN_CICE%PRES_DOCN%DOM_RTM_SGLC_SWAV -mach yellowstone -res f19_f19) with published number concentrations from Gettelman et al. (2008, http://acacia.ucar.edu/staff/andrew/papers/gettelman2008_micro.pdf). I'm getting much smaller values (see attached figures), and am wondering why that is.
Here's what I'm plotting:
- ignoring grid boxes with fractional cloud coverage below 0.0001 (or set it to 0.0001)
- divide numliq by cloud cover to get in-cloud number concentration, multiply by density to convert from 1/kg to 1/m3, divide by 1e6 to convert from 1/m3 to 1/cm3: number_concentration=numliq/cloud*rho/1e6
- zonal and monthly average (only over grid boxes with cloud cover above 0.0001)
- Spinup time is more than one year.
- If I do this online in a module or offline as post-processing yields the same result.
- For ice, the problem is similar (generally smaller by about a factor of 10)

Any comments would be greatly appreciated!
Thanks!
Andreas


Ferret code:
let cloud2=if (cloud gt 0.0001) then cloud
let concliq=numliq/cloud2*density
fill/vlimits=1000:0:100/lev=(40,180,10)(inf)/pal=inverse_grayscale concliq[x=@ave,l=@ave]/1e6
This should be equivalent to micro_mg_cam.F90:
icwnc(i,k) = state_loc%q(i,k,ixnumliq) / max(0.0001_r8,liqcldf(i,k)) * &
state_loc%pmid(i,k) / (287.15_r8*state_loc%t(i,k))

 

[santos]

It's hard to tell what's going on here without knowing a little more detail. Specifically, how do you get the cloud fraction, density, and numliq? If you are doing this in postprocessing, what history outputs are you using?My best guess is that you are using a different cloud fraction from the one in the code (the "liqcldf" in the code actually corresponds to the "AST" output).

 

[santos]

It's hard to tell what's going on here without knowing a little more detail. Specifically, how do you get the cloud fraction, density, and numliq? If you are doing this in postprocessing, what history outputs are you using?My best guess is that you are using a different cloud fraction from the one in the code (the "liqcldf" in the code actually corresponds to the "AST" output).

 

[andreasb]

Thanks for your comment! For the postprocessing, I was using the CLOUD, NUMLIQ, T (temperature), and pressure (lev) from the output, to calculate density rho=pressure[Pa]/(287.15*T), and now tried AST instead of CLOUD, which improves things a bit, but not enough, see atached new figure.Here's my online code:    call cnst_get_ind('NUMLIQ', ixnumliq)       cld_idx             = pbuf_get_index('CLD')     ! total cloud fraction    ast_idx            = pbuf_get_index('AST')

    call pbuf_get_field(pbuf, cld_idx,     cldn,   start=(/1,1,itim/), kount=(/pcols,pver,1/))
       rho(1:ncol,k)= &
              state%pmid(1:ncol,k)/(287.15_r8*state%t(1:ncol,k)   numliq_v(i,k)= state%q(i,k,ixnumliq) / max(0.0001_r8,cldn(i,k)) *rho(i,k)with cldn = AST or CLD thanks,Andreas

 

[andreasb]

Thanks for your comment! For the postprocessing, I was using the CLOUD, NUMLIQ, T (temperature), and pressure (lev) from the output, to calculate density rho=pressure[Pa]/(287.15*T), and now tried AST instead of CLOUD, which improves things a bit, but not enough, see atached new figure.Here's my online code:    call cnst_get_ind('NUMLIQ', ixnumliq)       cld_idx             = pbuf_get_index('CLD')     ! total cloud fraction    ast_idx            = pbuf_get_index('AST')

    call pbuf_get_field(pbuf, cld_idx,     cldn,   start=(/1,1,itim/), kount=(/pcols,pver,1/))
       rho(1:ncol,k)= &
              state%pmid(1:ncol,k)/(287.15_r8*state%t(1:ncol,k)   numliq_v(i,k)= state%q(i,k,ixnumliq) / max(0.0001_r8,cldn(i,k)) *rho(i,k)with cldn = AST or CLD thanks,Andreas

 

[andrew]

Hi Andreas, Thanks for the message. Sorry you are having trouble with the comparisons. The estimates from the paper were done with specific variables in CEM: AWNC/FREQL. This takes the number concentration with a threshold applied, averages it, and uses FREQL as a counter to determine how many points contain valid information. I also am not sure that the pressure correction should be applied to turn values into volume. The assumption is per meter cubed (m-3), so to get to cm-3, divide by 1.e^6: do not also assume that values are in kg units. See if that fixes the issue. We probably need to clarify units, but originally I tihnk they were desired as mixing ratio, but taht is not what the microphysics produces.  Cheers, Andrew. 

 

[andrew]

Hi Andreas, Thanks for the message. Sorry you are having trouble with the comparisons. The estimates from the paper were done with specific variables in CEM: AWNC/FREQL. This takes the number concentration with a threshold applied, averages it, and uses FREQL as a counter to determine how many points contain valid information. I also am not sure that the pressure correction should be applied to turn values into volume. The assumption is per meter cubed (m-3), so to get to cm-3, divide by 1.e^6: do not also assume that values are in kg units. See if that fixes the issue. We probably need to clarify units, but originally I tihnk they were desired as mixing ratio, but taht is not what the microphysics produces.  Cheers, Andrew. 

 

[santos]

Hmm, I think that multiplying by the density actually is correct. There are certain steps in the microphysics that work with number per cm^3, but they do so by converting on the fly; ncic and other such variables are per kg. In fact they really must be converted to concentrations per kg of air once they are put into the CAM state object, which is where Andreas is retrieving it from.Using AWNC/FREQL might be an interesting check, but at first glance it looks to me like it's equivalent.I wonder if there is an issue with the runs themselves not being comparable, i.e. there may be some difference due to the different parameterizations in CAM5, or different tuning parameters.

 

[santos]

Hmm, I think that multiplying by the density actually is correct. There are certain steps in the microphysics that work with number per cm^3, but they do so by converting on the fly; ncic and other such variables are per kg. In fact they really must be converted to concentrations per kg of air once they are put into the CAM state object, which is where Andreas is retrieving it from.Using AWNC/FREQL might be an interesting check, but at first glance it looks to me like it's equivalent.I wonder if there is an issue with the runs themselves not being comparable, i.e. there may be some difference due to the different parameterizations in CAM5, or different tuning parameters.

 

[andrew]

I checked, and AWNC is the in-cloud number concentration in #/m-3 (with density applied).  I would try using AWNC/FREQL and see if that gets the same answer: I believe that probably what is happening is you are including a lot of points with small amounts of liquid water: the microphysics does not clean all of these up, but AWNC has thresholds:liqcldf_grid(i,k) > 0.01_r8 .and. icwmrst_grid(i,k) > 5.e-5_r8


I suspect that you are including a lot of points with low cloud fraction or (especially) small water content (empty clouds) that then have small number.

You are using instantaneous (single time step, not averaged) output yes? With AWNC/FREQL you do not need to.

Cheers,

Andrew

 

[andrew]

I checked, and AWNC is the in-cloud number concentration in #/m-3 (with density applied).  I would try using AWNC/FREQL and see if that gets the same answer: I believe that probably what is happening is you are including a lot of points with small amounts of liquid water: the microphysics does not clean all of these up, but AWNC has thresholds:liqcldf_grid(i,k) > 0.01_r8 .and. icwmrst_grid(i,k) > 5.e-5_r8


I suspect that you are including a lot of points with low cloud fraction or (especially) small water content (empty clouds) that then have small number.

You are using instantaneous (single time step, not averaged) output yes? With AWNC/FREQL you do not need to.

Cheers,

Andrew

 

[andreasb]

Thanks Andrew and Sean!So in the output I find that AWNC = density * NUMLIQ as I think it should be. However, using FREQL instead of AST makes a huge difference, see attached for new number concentration plot, which is now pretty close to Andrew's original one! But look at the second plot - FREQL (liqcldf) and AST seem to be something very different, despite that it says in the code   ! Microphysics assumes 'liquid stratus frac = ice stratus frac
   !                      = max( liquid stratus frac, ice stratus frac )'.
   alst_mic => ast
   aist_mic => ast

and: if ( liqcldf(i,k) > 0.01_r8 .and. icwmrst(i,k) > 5.e-5_r8 ) then
  liqcldf(:ncol,top_lev:pver) = ast(:ncol,top_lev:pver)
... I also tested if FREQL+FREQI = AST (see attached), but that's not the case either. So could you elucidate me on what these actually mean? That would be very much appreciated.Thanks again for your help!Andreas  

 

[andreasb]

Thanks Andrew and Sean!So in the output I find that AWNC = density * NUMLIQ as I think it should be. However, using FREQL instead of AST makes a huge difference, see attached for new number concentration plot, which is now pretty close to Andrew's original one! But look at the second plot - FREQL (liqcldf) and AST seem to be something very different, despite that it says in the code   ! Microphysics assumes 'liquid stratus frac = ice stratus frac
   !                      = max( liquid stratus frac, ice stratus frac )'.
   alst_mic => ast
   aist_mic => ast

and: if ( liqcldf(i,k) > 0.01_r8 .and. icwmrst(i,k) > 5.e-5_r8 ) then
  liqcldf(:ncol,top_lev:pver) = ast(:ncol,top_lev:pver)
... I also tested if FREQL+FREQI = AST (see attached), but that's not the case either. So could you elucidate me on what these actually mean? That would be very much appreciated.Thanks again for your help!Andreas  

 

[andrew]

Hi Andreas,FREQL is NOT the cloud fraction. FREQL is the frequency of occurrence of liquid water above a certain threshold (AST > 0.01 and In-cloud mixing ratio > 5e.-5 kg/kg (50ppm)). Ice works the same way. These are variables which are set to 1 every time AWNC has a value, so when averaged over time (a month): AWNC/FREQL recovers the in cloud values only when clouds are present.Does that make sense? See here for some more explanation on the diagnostic fields:  http://www.cgd.ucar.edu/staff/andrew/diag_fields2.txtCheers, Andrew

 

[andrew]

Hi Andreas,FREQL is NOT the cloud fraction. FREQL is the frequency of occurrence of liquid water above a certain threshold (AST > 0.01 and In-cloud mixing ratio > 5e.-5 kg/kg (50ppm)). Ice works the same way. These are variables which are set to 1 every time AWNC has a value, so when averaged over time (a month): AWNC/FREQL recovers the in cloud values only when clouds are present.Does that make sense? See here for some more explanation on the diagnostic fields:  http://www.cgd.ucar.edu/staff/andrew/diag_fields2.txtCheers, Andrew

 

[andreasb]

Thanks Andrew, that helped a lot!Just to confirm that I use this correctly then, I would very much appreciate if you could tell me if I'm doing this correctly then. For calculating electrical conductivity in clouds, I require the spectral distribution of droplets and ice inside the cloud, using the gamma function as it's used in the code and described in Morrison and Gettelman 2008, eq. 1-4. See code below.Thanks a lot!Andreas              ! liquid:             Ndash=state%q(i,k,ixnumliq)/freql(i,k)
             qdash=state%q(i,k,ixcldliq)/freql(i,k)
             mu_l(i,k)=0.0005714_r8*(Ndash/1.e6_r8*rho(i,k))+0.2714_r8
             mu_l(i,k)=1._r8/(mu_l(i,k)**2)-1._r8
             mu_l(i,k)=max(mu_l(i,k),2._r8)
             mu_l(i,k)=min(mu_l(i,k),15._r8)
             lambda_l(i,k) = (pi/6._r8*rhow*Ndash*gamma(mu_l(i,k)+4._r8)/ &
                     (qdash*gamma(mu_l(i,k)+1._r8)))**(1._r8/3._r8)
             N0_l(i,k)=(Ndash*lambda_l(i,k)**(mu_l(i,k)+1.))/gamma(mu_l(i,k)+1._r8)
             ! ice
             Ndash=state%q(i,k,ixnumice)/freqi(i,k)
             qdash=state%q(i,k,ixcldice)/freqi(i,k)
             lambda_i(i,k) = (pi*rhoi*Ndash/qdash)**(1._r8/3._r8)
             N0_i(i,k)=Ndash*lambda_i(i,k)*2._r8
...    ! spectrum
    do j=1,numbins
       r_l(j) = j* 0.5e-6 !  up to 100 microns
       r_i(j) = j*1.25e-6 !  up to 250 microns
    end do! droplet size spectral distribution [#/m3/m]             nsl=rho(i,k)* N0_l(i,k)*(r_l(j)*2)**mu_l(i,k)*exp(-lambda_l(i,k)*r_l(j)*2)
             nsi=rho(i,k)* N0_i(i,k)                            *exp(-lambda_i(i,k)*r_i(j)*2)           

 

[andreasb]

Thanks Andrew, that helped a lot!Just to confirm that I use this correctly then, I would very much appreciate if you could tell me if I'm doing this correctly then. For calculating electrical conductivity in clouds, I require the spectral distribution of droplets and ice inside the cloud, using the gamma function as it's used in the code and described in Morrison and Gettelman 2008, eq. 1-4. See code below.Thanks a lot!Andreas              ! liquid:             Ndash=state%q(i,k,ixnumliq)/freql(i,k)
             qdash=state%q(i,k,ixcldliq)/freql(i,k)
             mu_l(i,k)=0.0005714_r8*(Ndash/1.e6_r8*rho(i,k))+0.2714_r8
             mu_l(i,k)=1._r8/(mu_l(i,k)**2)-1._r8
             mu_l(i,k)=max(mu_l(i,k),2._r8)
             mu_l(i,k)=min(mu_l(i,k),15._r8)
             lambda_l(i,k) = (pi/6._r8*rhow*Ndash*gamma(mu_l(i,k)+4._r8)/ &
                     (qdash*gamma(mu_l(i,k)+1._r8)))**(1._r8/3._r8)
             N0_l(i,k)=(Ndash*lambda_l(i,k)**(mu_l(i,k)+1.))/gamma(mu_l(i,k)+1._r8)
             ! ice
             Ndash=state%q(i,k,ixnumice)/freqi(i,k)
             qdash=state%q(i,k,ixcldice)/freqi(i,k)
             lambda_i(i,k) = (pi*rhoi*Ndash/qdash)**(1._r8/3._r8)
             N0_i(i,k)=Ndash*lambda_i(i,k)*2._r8
...    ! spectrum
    do j=1,numbins
       r_l(j) = j* 0.5e-6 !  up to 100 microns
       r_i(j) = j*1.25e-6 !  up to 250 microns
    end do! droplet size spectral distribution [#/m3/m]             nsl=rho(i,k)* N0_l(i,k)*(r_l(j)*2)**mu_l(i,k)*exp(-lambda_l(i,k)*r_l(j)*2)
             nsi=rho(i,k)* N0_i(i,k)                            *exp(-lambda_i(i,k)*r_i(j)*2)           

 

[andrew]

Hi Andreas, I cannot tell if this is on-line or off-line, but if it is off line you should NOT be using freqi and freql. They are just 0 or 1 in the code, and are only for output. I think you intend cloud fraction here. Andrew

 

[andrew]

Hi Andreas, I cannot tell if this is on-line or off-line, but if it is off line you should NOT be using freqi and freql. They are just 0 or 1 in the code, and are only for output. I think you intend cloud fraction here. Andrew

 

[i_h_h_karset@geo_uio_no]

Hi, Andrew. I can see that this link is old, but I'm also struggling with averaging CDNC (and the radius). I wonder if it's enough to divide by FREQL when finding the mean CDNC over time, lat, lev and lon, or if we also have to some how weight by the cloud fraction. My question is whether FREQL only contains information about how many of the timesteps the grid-box contained liquid (cloud), or if it also accounts for how much of the grid-box that is coved by a cloud. Thank you in advance for your help. RegardsInger Helene