Mersenne Twister in in-line post

[WenMeng-NOAA]

Dusan found the mersenne twister module in UPP routine CALWXT_BOURG.f are not working properly:

[149] incalwxtbg, rn 3.3978019E+38 -3.3952875E+38 [145] in calwxtbg, jsta,jend= 33 64 im= 384 [145] in calwxtbg,me= 1 iseed= 302134645 [145] incalwxtbg, rn 3.3978019E+38 -3.3952875E+38

The expected rn range should be within in 0 to 1. The mersenne_twister in UPP is from NCEPLIBS w3em. This issue is cased different library precision used in ufs forecast component ( w3emc_d) and upp lib (w3emc_4) for in-line post.

[WenMeng-NOAA]

Solution option 1) Dusan proposed the fix via replace subroutine random_number (from mersenne_twister module) with the one from standard fortran function instead of w3emc module. Dusan's proposal can be found at https://github.com/DusanJovic-NOAA/EMC_post/commit/b6aeef5cd12072827261d954b4eefee6d344b632

[WenMeng-NOAA]

Solution option 2) Moorthi indicated he had committed the fix for this kind of issue in his feature branch: https://github.com/SMoorthi-emc/EMC_post/commit/b87c8e8b91b19a69442d8f8be4642dd26bf91bbd

[WenMeng-NOAA]

Wen tested both option 1 and 2 with in-line post. Both option2 would fix rn value as expected range (1, 0). With runtime logs, the expected changes would be like from incalwxtbg, rn 3.3978019E+38 -3.3952875E+38 into incalwxtbg, rn 0.999998649116550 2.458319165258258E-005

[WenMeng-NOAA]

Huiya's comments:

I believe we switched UPP to use mersenne twister many years ago because the Fortran random number generator had a bias. This is consistent with what Jun said that the former is a better random number generator.

We're hoping to replace precipitation type algorithms with GSL one so calwxt_* will be potentially retired in GFS v17. With that said, I'm not sure how much time we want to spend on this.

[WenMeng-NOAA]

Wen created a new feature branch post_random_number and incorporate Moorthi's fixes at https://github.com/WenMeng-NOAA/EMC_post/commit/2327206d2db9558d0b5a868f914e8559d04a064a

[climbfuji]

My suggested approach for the mersenne twister is summarized in an issue for the ufs-weather-model: https://github.com/ufs-community/ufs-weather-model/issues/324

[junwang-noaa]

Wen, in general to make POST not fail at this time before the final mersenne twister interface is ready, I think you need to make sure POST gets correct random number for both real4 and real8 mersenne twister interfaces by checking the data range. So it is good to uncomment the changes you made in CALWXT_BOURG.f or make a small subroutine to handle this random number interface check.

On Thu, Dec 17, 2020 at 10:08 AM Dom Heinzeller notifications@github.com wrote:

My suggested approach for the mersenne twister is summarized in an issue for the ufs-weather-model: ufs-community/ufs-weather-model#324 https://github.com/ufs-community/ufs-weather-model/issues/324

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/EMC_post/issues/242#issuecomment-747495122, or unsubscribe https://github.com/notifications/unsubscribe-auth/AI7D6TNVZISAFXN2PQGNB4TSVINFZANCNFSM4U6OR7EQ .

[WenMeng-NOAA]

@junwang-noaa Yes, the commented block in CALWXT_BOURG.f as you suggested from our last discussion. I will uncomment that block for testing. There is related fixes in SURFACE.f. Do you have more comments on that part?

[junwang-noaa]

I think you can create a subroutine so that it can be called wherever mersenne twister random number is called. Basically the subroutine will return a correct real 4 random number for post when either real4 or real8 mersenne twister interface is used.

On Thu, Dec 17, 2020 at 10:48 AM WenMeng-NOAA notifications@github.com wrote:

@junwang-noaa https://github.com/junwang-noaa Yes, the commented block in CALWXT_BOURG.f as you suggested from our last discussion. I will uncomment that block for testing. There is related fixes in SURFACE.f. Do you have more comments on that part?

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

[WenMeng-NOAA]

Are you meaning a new interface in CALWXT_BOURG.f like

subroutine calwxt_bourg_post(im,jm,jsta_2l,jend_2u,jsta,jend,lm,lp1, & & iseed,g,pthresh, & & t,q,pmid,pint,lmh,prec,zint,ptype,me) ... real rn(imjm2) ... ... interface right_random_nuber subroutine random_number_4(rn) implicit none real(kind=4) :: rn(imjm2) call random_number(rn) end subroutine random_number_4 subroutine random_number_8(rn) implicit none real(kind=8) :: rn(imjm2) call random_number(rn) end subroutine random_number_8 end interface call right_random_num(rn) ... ... end subroutine calwxt_bourg_post

On Thu, Dec 17, 2020 at 10:56 AM Jun Wang notifications@github.com wrote:

I think you can create a subroutine so that it can be called wherever mersenne twister random number is called. Basically the subroutine will return a correct real 4 random number for post when either real4 or real8 mersenne twister interface is used.

On Thu, Dec 17, 2020 at 10:48 AM WenMeng-NOAA notifications@github.com wrote:

@junwang-noaa https://github.com/junwang-noaa Yes, the commented block in CALWXT_BOURG.f as you suggested from our last discussion. I will uncomment that block for testing. There is related fixes in SURFACE.f. Do you have more comments on that part?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://github.com/NOAA-EMC/EMC_post/issues/242#issuecomment-747523251 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AI7D6TJWSU6O4QNRR7LEMPTSVIR5ZANCNFSM4U6OR7EQ

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/EMC_post/issues/242#issuecomment-747528267, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALQGNEQC4VDEBYCIYVNSUBDSVISY5ANCNFSM4U6OR7EQ .

-- Wen Meng IMSG at NOAA/NWS/NCEP/EMC 5830 University Research Ct., Room 2070 College Park, MD 20740 Wen.Meng@noaa.gov 301-683-3779

[SMoorthi-emc]

Wen, It is possible that my update may not work in an offline post if mersenne-twister is compiled in single precision. Moorthi

On Thu, Dec 17, 2020 at 10:07 AM WenMeng-NOAA notifications@github.com wrote:

Wen created a new feature branch post_random_number and incorporate Moorthi's fixes at WenMeng-NOAA@2327206 https://github.com/WenMeng-NOAA/EMC_post/commit/2327206d2db9558d0b5a868f914e8559d04a064a

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/EMC_post/issues/242#issuecomment-747494425, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALLVRYRGQ5ROL45KDZKO4QTSVINCBANCNFSM4U6OR7EQ .

-- Dr. Shrinivas Moorthi Research Meteorologist Modeling and Data Assimilation Branch Environmental Modeling Center / National Centers for Environmental Prediction 5830 University Research Court - (W/NP23), College Park MD 20740 USA Tel: (301)683-3718

e-mail: Shrinivas.Moorthi@noaa.gov Phone: (301) 683-3718 Fax: (301) 683-3718

[WenMeng-NOAA]

Just commit https://github.com/WenMeng-NOAA/EMC_post/commit/585c91deb78ce7b0a32192600713d92f31f60c6d per Jun's last comment. A new interface for handle either real 4 or real 8. I am not familiar to interface structure in fortran. Please let me know your comments.

[junwang-noaa]

Wen, if I remember correctly, post can only run with real(4), so you don't need extra interface there. In generally if mersenne twister should provide such interface for the applications that calls it. SInce it is currently not available, what you need in the subroutine is something like: subroutine right_random_number(rn4) use mersenne_twister random_number real(kind=4), intent(out) :: rn4(imjm2)

real(kind=8) :: rn8(imjm2)

call random_number(rn4) if(rn4 > 1. or rn4<0.) then call random_number(rn8) rn4=rn8 if(rn4 > 1. or rn4<0.) print "error" endif end subroutine right_random_number

then you can call this right_random_numberin post CALWXT_BOURG.f or SURFACE.f. This could work for both standalone and inline post.

On Thu, Dec 17, 2020 at 1:28 PM WenMeng-NOAA notifications@github.com wrote:

Just commit WenMeng-NOAA@585c91d https://github.com/WenMeng-NOAA/EMC_post/commit/585c91deb78ce7b0a32192600713d92f31f60c6d per Jun's last comment. A new interface for handle either real 4 or real 8. I am not familiar to interface structure in fortran. Please let me know your comments.

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

[WenMeng-NOAA]

Jun,

Somehow, the job is crashed inside of the new subroutine right_random_number I add it in CALWXT_BOURG.f as:

  subroutine right_random_number(rn4,nm)
    use mersenne_twister, only: random_number
    implicit none
    integer,intent(in) :: nm
    real (kind=4),intent(out) :: rn4(nm)
    real (kind=8) :: rn8(nm)
    call random_number(rn4)
    if(maxval(rn4) > 1. .or. minval(rn4) < 0.) then
      call random_number(rn8)
      rn4=rn8
      if(maxval(rn4) > 1. .or. minval(rn4) < 0.) then
        print*, "in right_random_number, error"
      endif
    endif
   end subroutine right_random_number

---

My working version is at /scratch2/NCEPDEV/ovp/Wen.Meng/LIBS/src/upp_v10.0.1/EMC_post/sorc/ncep_post.fd. Do you think anything is missing?

Thanks,

Wen

On Thu, Dec 17, 2020 at 1:42 PM Jun Wang notifications@github.com wrote:

Wen, if I remember correctly, post can only run with real(4), so you don't need extra interface there. In generally if mersenne twister should provide such interface for the applications that calls it. SInce it is currently not available, what you need in the subroutine is something like: subroutine right_random_number(rn4) use mersenne_twister random_number real(kind=4), intent(out) :: rn4(imjm2)

real(kind=8) :: rn8(imjm2)

call random_number(rn4) if(rn4 > 1. or rn4<0.) then call random_number(rn8) rn4=rn8 if(rn4 > 1. or rn4<0.) print "error" endif end subroutine right_random_number

then you can call this right_random_numberin post CALWXT_BOURG.f or SURFACE.f. This could work for both standalone and inline post.

On Thu, Dec 17, 2020 at 1:28 PM WenMeng-NOAA notifications@github.com wrote:

Just commit WenMeng-NOAA@585c91d < https://github.com/WenMeng-NOAA/EMC_post/commit/585c91deb78ce7b0a32192600713d92f31f60c6d

per Jun's last comment. A new interface for handle either real 4 or real 8. I am not familiar to interface structure in fortran. Please let me know your comments.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://github.com/NOAA-EMC/EMC_post/issues/242#issuecomment-747617470 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AI7D6TMFFE3Z5FC44ZVXM3LSVJEVTANCNFSM4U6OR7EQ

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/EMC_post/issues/242#issuecomment-747624777, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALQGNETSN2F5UZHIREOS6ILSVJGI5ANCNFSM4U6OR7EQ .

-- Wen Meng IMSG at NOAA/NWS/NCEP/EMC 5830 University Research Ct., Room 2070 College Park, MD 20740 Wen.Meng@noaa.gov 301-683-3779

[WenMeng-NOAA]

The solution comes to add a new UPP module upp_right_mersenne_twister which includes random_number routines for real 4 and real 8 and two types of routines switch via the interface inside of the module. This solution was tested with in-line post, off-line post on Dell, Cray, Hera and Orion.

[WenMeng-NOAA]

I spent a lot of time to tweak Moorth's fixes (keeping mersenne twister from w3emc) for working in in-line and off-line post tests last week. From Dusan's comments on PR #243, it seems that random_number call in new module upp_right_mersenne_twister is from Fortran intrinsic subroutine other than Mersenne twister from w3emc (I don't have enough fortran engine knowledge for that kind of modification).

Dusan's fixes would be straightforward for my UPP tests. One concern is switching to generic fortran routine random_number might cause some scientific issue. That's the reason mersenne_twister was adapted in UPP code many years ago. As Huiya's UPP project plan, we would adapt GSD precipitation type calculation method in GLOBAL model process in the future.

I would wait for comments from @HuiyaChuang-NOAA , @junwang-noaa and @SMoorthi-emc which fixes option we should adapt.

[HuiyaChuang-NOAA]

I support Wen's solution to use Dusan's fix.

Please note that this fix will not be used in GFS V16 implementation and by the time RRFS is implemented, UPP will most likely already transition to GSL's explicit precipitation type algorithm. Once this fix is committed, it will support Public release.

Again, the subroutines calwxt* will become obsolete about a year from now. Wen already spent a lot of time on this. Let's use Dusan's solution and move on.

Huiya