-
Notifications
You must be signed in to change notification settings - Fork 9
Remove inconsistencies in computation of air density with Thompson MP #79
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Changes from 6 commits
5f7968b
a0b31a5
9cd399d
a2d800b
5a9d256
8449be5
4c74498
f06e592
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -147,11 +147,11 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do | |
| ! Cloud particle sizes and number concentrations... | ||
|
|
||
| ! First, prepare cloud mixing-ratios and number concentrations for Calc_Re | ||
| rho = p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev)) | ||
| orho = 1./rho | ||
| do iLay = 1, nLev | ||
| do iCol = 1, nCol | ||
| qv_mp(iCol,iLay) = q_lay(iCol,iLay)/(1.-q_lay(iCol,iLay)) | ||
| rho(iCol,iLay) = 0.622*p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev)*(qv_mp(iCol,iLay)+0.622)) | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. To save some memory, shouldn't we be using scalar (singular) values inside the loops not arrays? Unless we are saving them for later usage, in which case I guess it's necessary. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think the loops in lines 149-166 and in lines 169-178 can be combined, which would allow us to make rho and orho scalars.
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I added here just one line at request from Dustin Swales. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Note that this line is also incorrect, because it multiplies arrays
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Of course! A copy/paste problem Thank you for fixing this. |
||
| orho(iCol,iLay) = 1./rho(iCol,iLay) | ||
| qc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq) / (1.-q_lay(iCol,iLay)) | ||
| qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay)) | ||
| qs_mp(iCol,iLay) = tracer(iCol,iLay,i_cldsnow) / (1.-q_lay(iCol,iLay)) | ||
|
|
@@ -169,7 +169,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do | |
| do iLay = 1, nLev | ||
| do iCol = 1, nCol | ||
| if (ltaerosol .and. qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then | ||
| nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho(iCol,iLay), nwfa(iCol,iLay)) * orho(iCol,iLay) | ||
| nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho(iCol,iLay), nwfa(iCol,iLay)*rho(iCol,iLay)) * orho(iCol,iLay) | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Same comment as before. This becomes a 3rd time for possible differences in code to exist. I don't see it as necessary. Is this possibly because timestep#1 itself calls radiation before microphysics? If so, we should resolve this problem permanently like it can be done in WRF using module_initialize_real.F to ensure no mass can exist without a corresponding number for any MP scheme for any species in 2 moments. It would be far more general and safe that way. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I don't recall the reasoning, but it shouldn't be the timestep, because There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. To address your other comment about 3rd time. If it isn't necessary for rrtmgp, then it isn't necessary for rrtmg and will disappear in both. Let's hear what the GSL developers have to say. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I don't know the history why this call to make_dropletnumber() is in here. My assumption was that the #-concentration needed to be updated consistently before calling radiation?
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I commented earlier on make_DropletNumber. We have to do it every time step for subgrid clouds.
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @dustinswales Your assumption is correct.
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The rrtmg_pre code prepares effective radii for the calls to progcld* for all microphysics schemes. Therefore it is logical to call calc_effecRad for THompson MP in this subroutine. RRTMGP has a separate subroutine for Thompson MP: GFS_rrtmgp_thompsonmp_pre.F90. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Hmm, you are right. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I strongly urge the group to consider coordinating a single place where we treat creation of water droplet number or ice crystal number for sub-grid-scale clouds in a unified single location, not among one or more variants of radiation physics. I am already handling explicit/grid-scale clouds within my MP scheme in a uniform way. We should do likewise for SGS clouds. And, in so doing, we can also permit potential other SGS cloud schemes to do likewise. In my view, this should be really simple to do in the file Dustin mentioned (module_SGSCloud_RadPre.F90). Could you also look over my code mods in PR#567 and note that I have included specific size constraints as parameter constants at top of MP-thompson that could be linked via USE ..., ONLY .... statements? And, in a related note, my changes to the radiative effective radii also removed the internal size constraints such that after calling the size calculation, the MIN/MAX statements are now quite important. This permits the actual calculation within the subroutine to go beyond the lower/upper limits in contrast to the old behavior. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @gthompsnJCSDA Conceptually, the SGS clouds are not related to this step in the process. The SGS clouds come in from other physical parameterizations (e.g. BL and convection) and the role of the sgsCloud_radpre routine is to couple these other MP to the cloud fields, which needs to be done before we couple the MP to the radiation. I see no problem with importing min/max's if there needs to be bounding on the effective radii. I do exactly that for the RRTMGP particle sizes. |
||
| endif | ||
| if (qi_mp(iCol,iLay) > 1.e-12 .and. ni_mp(iCol,iLay) < 100.) then | ||
| ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho(iCol,iLay), t_lay(iCol,iLay)) * orho(iCol,iLay) | ||
|
|
||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -692,7 +692,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to | |
| ! local variables | ||
| integer :: i,k,n,tracers | ||
|
|
||
| real(kind=kind_phys), dimension(im,levs) :: rho_dryair | ||
| real(kind=kind_phys), dimension(im,levs) :: rho_air | ||
| real(kind=kind_phys), dimension(im,levs) :: qv_mp !< kg kg-1 (dry mixing ratio) | ||
| real(kind=kind_phys), dimension(im,levs) :: qc_mp !< kg kg-1 (dry mixing ratio) | ||
| real(kind=kind_phys), dimension(im,levs) :: qi_mp !< kg kg-1 (dry mixing ratio) | ||
|
|
@@ -742,16 +742,16 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to | |
| if (imp_physics == imp_physics_thompson .and. (ntlnc>0 .or. ntinc>0)) then | ||
| do k=1,levs | ||
| do i=1,im | ||
| !> - Density of air in kg m-3 | ||
| rho_dryair(i,k) = prsl(i,k) / (con_rd*save_tcp(i,k)) | ||
| !> - Convert specific humidity to dry mixing ratio | ||
| qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) | ||
| !> - Density of air in kg m-3 | ||
| rho_air(i,k) = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Same comment about memory footprint and possibly using a local scalar variable not arrays. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes, in this place certainly possible. One could also define another scalar orho = 1/rho and do only one division instead of two (754, 763). |
||
| if (ntlnc>0) then | ||
| !> - Convert moist mixing ratio to dry mixing ratio | ||
| qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) | ||
| !> - Convert number concentration from moist to dry | ||
| nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) | ||
| nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho_dryair(i,k), nwfa(i,k)) * (one/rho_dryair(i,k))) | ||
| nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho_air(i,k), nwfa(i,k)*rho_air(i,k)) * (one/rho_air(i,k))) | ||
| !> - Convert number concentrations from dry to moist | ||
| gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) | ||
| endif | ||
|
|
@@ -760,7 +760,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to | |
| qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) | ||
| !> - Convert number concentration from moist to dry | ||
| ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) | ||
| ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho_dryair(i,k), save_tcp(i,k)) * (one/rho_dryair(i,k))) | ||
| ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho_air(i,k), save_tcp(i,k)) * (one/rho_air(i,k))) | ||
| !> - Convert number concentrations from dry to moist | ||
| gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) | ||
| endif | ||
|
|
||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -159,10 +159,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & | |
| ! Geopotential height in m2 s-2 to height in m | ||
| hgt = phil/con_g | ||
|
|
||
| ! Density of air in kg m-3 and inverse density of air | ||
| rho = prsl/(con_rd*tgrs) | ||
| orho = 1.0/rho | ||
|
|
||
| ! Prior to calling the functions: make_DropletNumber, make_IceNumber, make_RainNumber, | ||
| ! the incoming mixing ratios should be converted to units of mass/num per cubic meter | ||
| ! rather than per kg of air. So, to pass back to the model state variables, | ||
|
|
@@ -178,6 +174,10 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & | |
| qs_mp = qs/(1.0_kind_phys-spechum) | ||
| qg_mp = qg/(1.0_kind_phys-spechum) | ||
|
|
||
| ! Density of air in kg m-3 and inverse density of air | ||
| rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) | ||
| orho = 1.0/rho | ||
|
|
||
| !> - Convert number concentrations from moist to dry | ||
| ni_mp = ni/(1.0_kind_phys-spechum) | ||
| nr_mp = nr/(1.0_kind_phys-spechum) | ||
|
|
@@ -304,7 +304,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & | |
|
|
||
| ! If qc is in boundary conditions but nc is not, calculate nc from qc, rho and nwfa | ||
| if (maxval(qc_mp)>0.0 .and. maxval(nc_mp)==0.0) then | ||
| nc_mp = make_DropletNumber(qc_mp*rho, nwfa) * orho | ||
| nc_mp = make_DropletNumber(qc_mp*rho, nwfa*rho) * orho | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This should already be a part of PR#567 There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. If it is, I'll be giving a merge conflict to resolve, so no problem. |
||
| end if | ||
|
|
||
| ! If nc is in boundary conditions but qc is not, reset nc to zero | ||
|
|
@@ -428,6 +428,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & | |
|
|
||
| ! Air density | ||
| real(kind_phys) :: rho(1:ncol,1:nlev) !< kg m-3 | ||
| !rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) | ||
| ! Hydrometeors | ||
| real(kind_phys) :: qv_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) | ||
| real(kind_phys) :: qc_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) | ||
|
|
@@ -510,7 +511,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & | |
| end if | ||
|
|
||
| !> - Density of air in kg m-3 | ||
| rho = prsl/(con_rd*tgrs) | ||
| rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) | ||
|
|
||
| !> - Convert omega in Pa s-1 to vertical velocity w in m s-1 | ||
| w = -omega/(rho*con_g) | ||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There are inconsistencies with how I am testing the limits of QC and NC in mp_thompson code. And why are we needing to make droplet number in radiation code? Shouldn't this only be done in one place (MP code) and the resulting number just be given to radiation schemes?
I don't like all the duplication of code for fear that something goes awry/different in one than another place. If we need to make another isolation of ensuring species number and mass jointly make sense, then let's do it once and never more than once. But radiation code seems weird to me calling droplet or ice number fixes.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@joeolson42 @tanyasmirnova @hannahcbarnes was the reason for this the logic in the
sgscloud_radpreandsgscloud_radpostschemes, i.e. the modification ofqcandqibefore radiation with boundary layer clouds from MYNN PBL insgscloud_radpre, and then restoring the old values in ``sgscloud_radpost`?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We have to add the call to make_DropletNumber here only for the case when there are subgrid clouds created in PBL and convective schemes. They create sub-grid scale mixing ratios, but no NCs. After that the calc_effectRad is called using consistent mixing rations and nc_mp and ni_mp. Effective radii are needed for the radiation. After call to radiation, the sub-grid clouds are subtracted to restore the resolved mixing ratios from Thompson MP.