Skip to content
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

Transition to Thompson Microphysics for Microwave All-sky Assimilation #743

Draft
wants to merge 19 commits into
base: develop
Choose a base branch
from

Conversation

azadeh-gh
Copy link

@azadeh-gh azadeh-gh commented May 7, 2024

Description
Transition from using the GFDL microphysics scheme, currently in operation, to the Thompson scheme (GFSv17) is proposed for the microwave all-sky assimilation. This transition aims to enhance the accuracy and reliability of forecasting. Issue#719

Resolves #719

Type of change

Please delete options that are not relevant.

  • Bug fix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)
  • Breaking change (fix or feature that would cause existing functionality to not work as expected)
  • This change requires a documentation update

How Has This Been Tested?

Checklist

  • My code follows the style guidelines of this project
  • I have performed a self-review of my own code
  • I have commented my code, particularly in hard-to-understand areas
  • New and existing tests pass with my changes
  • Any dependent changes have been merged and published

Copy link
Contributor

@RussTreadon-NOAA RussTreadon-NOAA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Check for GSI Code Standard compliance. Minor items noted. Did not fully check crtm_interface.f90.

src/gsi/crtm_interface.f90 Outdated Show resolved Hide resolved
src/gsi/crtm_interface.f90 Outdated Show resolved Hide resolved
src/gsi/crtm_interface.f90 Outdated Show resolved Hide resolved
src/gsi/crtm_interface.f90 Outdated Show resolved Hide resolved
src/gsi/crtm_interface.f90 Outdated Show resolved Hide resolved
src/gsi/crtm_interface.f90 Outdated Show resolved Hide resolved
src/gsi/crtm_interface.f90 Outdated Show resolved Hide resolved
src/gsi/general_read_gfsatm.f90 Outdated Show resolved Hide resolved
src/gsi/netcdfgfs_io.f90 Outdated Show resolved Hide resolved
src/gsi/netcdfgfs_io.f90 Outdated Show resolved Hide resolved
Copy link
Contributor

@RussTreadon-NOAA RussTreadon-NOAA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you, @azadeh-gh , for updating the code. It looks much better now in terms of GSI code standards.

+ sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
+ sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
smoc = a_ * smob**b_
reff(k) = MAX(2.51E-6_r_kind, MIN(0.5*(smoc/smob), 1999.E-6_r_kind))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

0.5 should be 0.5_r_kind

@RussTreadon-NOAA
Copy link
Contributor

Rebuilt azadeh-gh:feature/ThompsonMP_Azadeh at c24e568 on Cactus and rerun ctests. Channels 1 to 13 of gmi_gpm now also shows contrl versus updat differences. AMSU-A differences occur in channels 1 to 5 and 15. ATMS differences occur in channels 1 to 6 and 16 to 22.

@@ -1730,6 +1761,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
ges_qsat (ixp,iy ,k, itsigp)*w10+ &
ges_qsat (ix ,iyp,k, itsigp)*w01+ &
ges_qsat (ixp,iyp,k, itsigp)*w11)*dtsigp
rh(k) = q(k)/qsat(k) !added for cloud fraction computation
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@azadeh-gh The variable qsat is declared but not defined.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think qsat should qs

@@ -1730,6 +1761,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
ges_qsat (ixp,iy ,k, itsigp)*w10+ &
ges_qsat (ix ,iyp,k, itsigp)*w01+ &
ges_qsat (ixp,iyp,k, itsigp)*w11)*dtsigp
rh(k) = q(k)/qsat(k) !added for cloud fraction computation
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think qsat should qs

cf_calc = zero
! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr)
qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5)
call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qsat,cf_calc)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here as above for qsat.

Comment on lines 2131 to 2139
if (trim(cloud_names_fwd(ii))=='ni' .and. atmosphere(1)%temperature(k)<t0c) then
cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii))
cloud_efr(k,ii)=max(5.001_r_kind, cloud_efr(k,ii))
endif
if (trim(cloud_names_fwd(ii))=='nr' .and. atmosphere(1)%temperature(k)<t0c) then
cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii))
cloud_efr(k,ii)=max(5.001_r_kind, cloud_efr(k,ii))
endif

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These lines are not needed and should be removed.
We do not need to check number concentration since CRTM does not take number concentration as input. We just need number concentration (ni and nr) for effective radius calculation

@emilyhcliu
Copy link
Contributor

@azadeh-gh So far, we have made sure the code can handle imp_physics=8 and imp_imp_physics=11.
The last thing to check is the effective radius and cloud fraction from Thompson and GFDL.
I ran both cases and printed out all atmospheric profiles (T, q, ql, qi, qr, qs, qg, cloud fraction and effective radius of all hydrometeors) from both methods (Thompson and GFDL). I do not think the effective radius from the Thompson scheme is calculated correctly. Please check the subroutines (calc_thompson_reff and calc_thompson_cloudfrac) again.

They can be found on HERA:
/scratch2/NCEPDEV/stmp1/Emily.Liu/GSI382/thompson-debug-amsua
/scratch2/NCEPDEV/stmp1/Emily.Liu/GSI382/gfdl-debug-amsua

Look at the 10th profile in fort.1004102 for both cases (search icmask in the file until you see the first icmask is T (true))

 icmask        :  T
 lat lon       :   -57.3952000000000        106.477000000000
 channel no.   :            2
 emissivity    :   0.527780801232015        206.827145270325
tsim tsim_clr :    1.747967983657567572E+02   1.614273477412949944E+02   8.695621635317450915E-01
 tcc           :   0.869562163531745

It should be very easy to find the bug since we have GFDL as a reference. The values (effective radius) from Thompson should not be too different.

@emilyhcliu
Copy link
Contributor

emilyhcliu commented May 31, 2024

@azadeh-gh I performed a sanity check for two branches: develop (control) and feature/ThompsonMP_Azadeh (experiment)

GSI Builds, Scripts, and IC for the Sanity Check
The set up can be found on HERA:
control: /scratch1/NCEPDEV/da/Emily.Liu/git/GSI/GSI-emilyhcliu/GSI
experiment: /scratch1/NCEPDEV/da/Emily.Liu/git/GSI/GSI-emilyhcliu/GSI-thompson-azadeh

The scripts to run the two experiments:
/scratch1/NCEPDEV/da/Emily.Liu/git/GSI/GSI-emilyhcliu/scripts/rungsi_develop.sh
/scratch1/NCEPDEV/da/Emily.Liu/git/GSI/GSI-emilyhcliu/scripts/rungsi_thompson_azadeh_for_develop.sh

The initial conditions:
COMROOT=/scratch2/NCEPDEV/stmp3/Azadeh.Gholoubi/TestExperiment/MonthThompsonMP

Sanity Check Configuration
The sanity check has the following configurations for 2023050200 and identical results are expected.

  1. same initial condition (with the Thompson forecast; so the cloud profile has number concentrations for rain (nr) and ice (ni))
  2. same set of observing system
  3. 4denvar run
  4. default global anavinfo for the control; update global anavinfo (with nr and ni added) for the experiment
  5. imp_physics is set to 11 (GFDL) for both control and experiment
  6. use GFDL cloud table for both control and experiment

Test Run and Output Directory
Control: /scratch2/NCEPDEV/stmp1/Emily.Liu/GSI382/develop
Experiment: /scratch2/NCEPDEV/stmp1/Emily.Liu/GSI382/thompson-azadeh-for-develop

Results
I checked the Jo tables and the minimization step output for all 202 steps. They are identical between the Control and Experiment (GSI develop vs. your branch). The evidence is attached below.
The regression test for your branch against the develop should be passed.

Control at the end of the first minimization loop/ second outer loop

 36:  Begin J table inner/outer loop           0           2
 36:     J contribution                 1             2             3             4             5             6             7
 36: background              1.008904E+05  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
 36: excess   moisture       1.369165E-01  1.374345E-01  1.363932E-01  1.343722E-01  1.326088E-01  1.305097E-01  1.280860E-01
 36: surface pressure        1.914876E+03  1.204426E+03  1.277209E+03  2.650365E+03  1.194631E+03  1.212828E+03  1.259412E+02
 36: temperature             2.554828E+03  5.320637E+03  8.836285E+03  1.734312E+04  1.034742E+04  5.817320E+03  3.132768E+03
 36: wind                    8.087878E+03  2.197073E+04  2.761711E+04  4.641993E+04  3.415848E+04  2.043609E+04  9.243730E+03
 36: moisture                1.995685E+02  5.340923E+02  1.496204E+03  2.176303E+03  4.961715E+02  3.997539E+02  4.190078E+02
 36: sst                     3.793319E+02  3.853462E+02  3.804480E+02  4.858384E+02  4.747631E+02  4.413564E+02  5.844765E+01
 36: ozone                   1.391325E+03  2.662759E+03  2.652272E+03  1.791333E+03  2.379520E+03  2.215390E+03  1.913433E+03
 36: gps bending angle       1.715324E+04  3.820486E+04  3.878015E+04  3.878843E+04  3.869788E+04  4.171309E+04  2.063999E+04
 36: radiance                8.382490E+04  1.490322E+05  1.658015E+05  1.494925E+05  1.490969E+05  1.598365E+05  7.196505E+04
 36:      J term                                     J
 36: background                   1.0089041321553968E+05
 36: excess   moisture            9.3632088894809418E-01
 36: surface pressure             9.5802762609576021E+03
 36: temperature                  5.3352382572843781E+04
 36: wind                         1.6793395224131059E+05
 36: moisture                     5.7211010589933849E+03
 36: sst                          2.6055316431712367E+03
 36: ozone                        1.5006030751918539E+04
 36: gps bending angle            2.3397761775828814E+05
 36: radiance                     9.2904948129338946E+05
 36:  -----------------------------------------------------
 36:  J Global                    1.5181177231173015E+06

Experiment at the end of the first minimization loop / second outerloop

4452  36:  Begin J table inner/outer loop           0           2
4453  36:     J contribution                 1             2             3             4             5             6             7
4454  36: background              1.008904E+05  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
4455  36: excess   moisture       1.369165E-01  1.374345E-01  1.363932E-01  1.343722E-01  1.326088E-01  1.305097E-01  1.280860E-01
4456  36: surface pressure        1.914876E+03  1.204426E+03  1.277209E+03  2.650365E+03  1.194631E+03  1.212828E+03  1.259412E+02
4457  36: temperature             2.554828E+03  5.320637E+03  8.836285E+03  1.734312E+04  1.034742E+04  5.817320E+03  3.132768E+03
4458  36: wind                    8.087878E+03  2.197073E+04  2.761711E+04  4.641993E+04  3.415848E+04  2.043609E+04  9.243730E+03
4459  36: moisture                1.995685E+02  5.340923E+02  1.496204E+03  2.176303E+03  4.961715E+02  3.997539E+02  4.190078E+02
4460  36: sst                     3.793319E+02  3.853462E+02  3.804480E+02  4.858384E+02  4.747631E+02  4.413564E+02  5.844765E+01
4461  36: ozone                   1.391325E+03  2.662759E+03  2.652272E+03  1.791333E+03  2.379520E+03  2.215390E+03  1.913433E+03
4462  36: gps bending angle       1.715324E+04  3.820486E+04  3.878015E+04  3.878843E+04  3.869788E+04  4.171309E+04  2.063999E+04
4463  36: radiance                8.382490E+04  1.490322E+05  1.658015E+05  1.494925E+05  1.490969E+05  1.598365E+05  7.196505E+04
4464  36:      J term                                     J
4465  36: background                   1.0089041321553968E+05
4466  36: excess   moisture            9.3632088894809418E-01
4467  36: surface pressure             9.5802762609576021E+03
4468  36: temperature                  5.3352382572843781E+04
4469  36: wind                         1.6793395224131059E+05
4470  36: moisture                     5.7211010589933849E+03
4471  36: sst                          2.6055316431712367E+03
4472  36: ozone                        1.5006030751918539E+04
4473  36: gps bending angle            2.3397761775828814E+05
4474  36: radiance                     9.2904948129338946E+05
4475  36:  -----------------------------------------------------
4476  36:  J Global                    1.5181177231173015E+06

Here are the final minimization step evaluation:
Control

cost,grad,step,b,step? =   2 150  1.465160993399310159E+06  1.971038208820891313E+01  2.822385197809975099E 01  1.604897083342762043E+00  good

Experiment

cost,grad,step,b,step? =   2 150  1.465160993399310159E+06  1.971038208820891313E+01  2.822385197809975099E-01  1.604897083342762043E+00  good

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Transition to Thompson Microphysics Scheme for Microwave All-sky Assimilation
3 participants