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
DusanJovic-NOAA@b6aeef5

[WenMeng-NOAA]

Solution option 2)
Moorthi indicated he had committed the fix for this kind of issue in his feature branch:
SMoorthi-emc@b87c8e8

[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
WenMeng-NOAA@2327206

[climbfuji]

My suggested approach for the mersenne twister is summarized in an issue for the ufs-weather-model: ufs-community/ufs-weather-model#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 ***@***.***> wrote: My suggested approach for the mersenne twister is summarized in an issue for the ufs-weather-model: ufs-community/ufs-weather-model#324 <ufs-community/ufs-weather-model#324> — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#242 (comment)>, 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 ***@***.***> 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 <#242 (comment)>, 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(im*jm*2) ... ... interface right_random_nuber subroutine random_number_4(rn) implicit none real(kind=4) :: rn(im*jm*2) call random_number(rn) end subroutine random_number_4 subroutine random_number_8(rn) implicit none real(kind=8) :: rn(im*jm*2) 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 ***@***.***> 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 ***@***.***> 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 > <#242 (comment) >, > 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 <#242 (comment)>, 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 ***@***.***> wrote: Wen created a new feature branch post_random_number and incorporate Moorthi's fixes at ***@***.*** <WenMeng-NOAA@2327206> — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#242 (comment)>, 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 WenMeng-NOAA@585c91d 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(im*jm*2) real(kind=8) :: rn8(im*jm*2) 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 ***@***.***> wrote: Just commit ***@***.*** <WenMeng-NOAA@585c91d> 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 <#242 (comment)>, 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 ***@***.***> 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(im*jm*2) real(kind=8) :: rn8(im*jm*2) 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 ***@***.***> wrote: > Just commit ***@***.*** > < WenMeng-NOAA@585c91d > > 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 > <#242 (comment) >, > 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 <#242 (comment)>, 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