# Weather and climate modeling codes from Fortran to C++ **URL:** \ **Category:** Uncategorized\ **Created:** [July 9, 2021, 2:30am UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531 "2021-07-09T02:30:55Z")\ **Posts on this page:** 20\ **Page:** 3
**Author:** ![fedebenelli](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/fedebenelli/32/1823_2.png) [@fedebenelli](https://fortran-lang.discourse.group/u/fedebenelli)\ **Post date:** [March 26, 2024, 2:24pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/42 "2024-03-26T14:24:07Z")
A more user-friendly interface for LAPACK/BLAS/etc libraries, I only use it for solving linear systems (for now at least) and when starting to learn with inherited codes I couldn’t find much sense to that “weird” `dgsev` call until reading further. Now my students have the wrapper function `solve_linear_system` So the call ```fortran call dgesv(n, nrhs, A, lda, ipiv, b, ldb, info) ``` becomes ```fortran X = solve_linear_system(A, b) ``` This is a really simple/basic example, but I think that correctly wrapping this routines (being more similar to numpy in some way) is a good approach to get a better Fortran ecosystem ---
**Author:** ![fedebenelli](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/fedebenelli/32/1823_2.png) [@fedebenelli](https://fortran-lang.discourse.group/u/fedebenelli)\ **Post date:** [March 26, 2024, 4:48pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/43 "2024-03-26T16:48:10Z")
I see that @FedericoPerini wrapper implementation is already doing this, makes me really happy! ---
**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\ **Post date:** [March 26, 2024, 5:25pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/44 "2024-03-26T17:25:17Z")
@fedebenelli that is exactly right. The high level design: - if you just want to solve a system (and stop the program if it can’t be solved), `x = solve(A, b)` will work - if you want to handle errors yourself, you can do `x = solve(A, b, err)` and then you can handle `err` yourself and have full flexibility. ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [March 26, 2024, 5:46pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/45 "2024-03-26T17:46:10Z")
I think the 2nd case should be written `call solve_system(a, b, x, err)` although it’s not a big deal, since one can easily write a subroutine wrapper to do this. Richard Maine (editor of the Fortran 2003 standard) [explained](https://comp.lang.fortran.narkive.com/AgDI7u8g/error-condition-testing) in comp.lang.fortran 16 years ago why he avoids having functions return error flags. > I’m even less a fan of having functions return values in arguments. > I’d probably recast the function as a subroutine if an error return was > a possibility. After all, to me most of the point of functions is to be > able to use them in the middle of expressions, and if your function is > instead invoked in the middle of an expression such as > > something = other\_stuff + whatever\*evaluate(xin) > > then you’ve “lost” your opportunity to test for the error. It is both > too late in that you have already used the “bad” result, plus the > appropriate test on something is now non-trivial. That’s probably a > bigger caveat than the others mentioned above, but I also presume it is > obvious. In my view, if a function can’t be used in an expression like > this, then there isn’t much point in it being a function. But that’s > more my personal style choice than a real caveat. ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [March 26, 2024, 5:57pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/46 "2024-03-26T17:57:13Z")
It will be impossible to make everybody happy, but we’ve been discussing interfaces at length with @jeremie.vandenplas, @everythingfunctional, @hkvzjal, @gnikit et al. during our past Monthly calls. So far we’ve come up with some interesting things: - All procedures that could be made `pure` were made `pure` (first time ever?) - Generalised `KIND`-dependent functions and 128-bit linear algebra support (first time ever?) - 2/3 interfaces for each operation, ideally: We also do not want to deviate too much from common interfaces such as NumPy and SciPy, because that’s what most people are used to nowadays. See for example the current matrix inversion interface: > Of course all comments and ideas are very welcome and we’ll do our best to optimize across all the requests! ---
**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\ **Post date:** [March 26, 2024, 6:44pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/47 "2024-03-26T18:44:19Z")
> [@Beliavsky](#): > > I think the 2nd case should be written > > `call solve_system(a, b, x, err)` > > although it’s not a big deal, since one can easily write a subroutine wrapper to do this. > > Richard Maine (editor of the Fortran 2003 standard) [explained](https://comp.lang.fortran.narkive.com/AgDI7u8g/error-condition-testing) in comp.lang.fortran 16 years ago why he avoids having functions return error flags. I fully agree with that. The other drawback of a return flags is that the function cannot be `pure` (or `simple` in F2023), which potentially prevent some compiler optimisations. In the current design choice purity is apparently reserved for the operator flavor, but not all operations can have an operator version. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [March 26, 2024, 6:50pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/48 "2024-03-26T18:50:14Z")
> [@FedericoPerini](#): > > Generalised `KIND`-dependent functions and 128-bit linear algebra support (first time ever?) Nice to hear about the progress you’re making with user-friendly linear algebra. Concerning high-precision LAPACK and BLAS, there is at least one implementation from a Japanese group([GitHub - nakatamaho/mplapack: The MPLAPACK: multiple precision version of BLAS and LAPACK](https://github.com/nakatamaho/mplapack)). I think I’ve seen other high-precision implementations before but can’t find the link now. ---
**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\ **Post date:** [March 26, 2024, 10:12pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/49 "2024-03-26T22:12:39Z")
It looks like @FedericoPerini has the subroutine version, as requested. However, it is not pure. Why not? ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [March 27, 2024, 7:24am UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/50 "2024-03-27T07:24:54Z")
Some of the high-level APIs have not been `pure`ified yet because we’re still focussing on the BLAS/LAPACK backends, but the will come soon! However, some functions in LAPACK contain `intent(out)` arguments, which prevents them, and all the upstream procedures using them, from being `pure`. So, some operations i.e. SVD cannot be made pure unless we change the LAPACK interface (which would break compatibility with the external libraries, or require large workarounds). One random example: > ---
**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\ **Post date:** [March 27, 2024, 9:06am UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/52 "2024-03-27T09:06:33Z")
> [@FedericoPerini](#): > > However, some functions in LAPACK contain `intent(out)` arguments, which prevents them, and all the upstream procedures using them, from being `pure`. Yes, but it’s a separate issue. Maybe in some future this will be fixed, by changing the interfaces of these functions (and of the underlying implementation). In contrast, there will be no way to make the higher level functions pure in the future if they have an intent(out) return flag. ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [March 27, 2024, 12:39pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/53 "2024-03-27T12:39:08Z")
> [@PierU](#): > > the higher level functions for those relying on `pure` BLAS/LAPACK procedures 2 out of 3 interfaces can be pure: - subroutine interface - `.operator.` interface ---
**Author:** ![rwmsu](https://avatars.discourse-cdn.com/v4/letter/r/48db29/32.png) [@rwmsu](https://fortran-lang.discourse.group/u/rwmsu)\ **Post date:** [March 27, 2024, 2:59pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/54 "2024-03-27T14:59:41Z")
I’ve hoped for a long time that the standards committee would allow programmers to use any ASCII symbol they wanted in a user defined operator so I could mimic MATLAB and use a backslash () as an operator to solve linear systems. ie ```auto interface operator(.\.) function solve (a, b) result(x) .. end function end interface x = a.\.b ``` I know it’s never going to happen but its something I would probably use a lot if available. ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [March 27, 2024, 3:16pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/55 "2024-03-27T15:16:11Z")
> [@rwmsu](#): > > ```auto > interface operator(.\.) > function solve (a, b) result(x) > .. > end function > end interface > > x = a.\.b > > ``` I’m not opposing your idea, but we could define a `.div.` operator and write x = `b.div.a` ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [March 27, 2024, 6:37pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/56 "2024-03-27T18:37:37Z")
Great idea! Let’s find one that would be suitable for LU solve in stdlib. Some further ideas: - `.backslash.` (I would rather avoid `.bs.` … lol) - `.ld.` (left divide? like Matlab’s mldivide) - `.inv.` that would support both `.inv.A` meaning A^{-1} and `A.inv.B` meaning A^{-1}B - `//` (yes, we may overload the intrinsic concatenation, looks close to a division but may be confusing?) ---
**Author:** ![rwmsu](https://avatars.discourse-cdn.com/v4/letter/r/48db29/32.png) [@rwmsu](https://fortran-lang.discourse.group/u/rwmsu)\ **Post date:** [March 27, 2024, 6:48pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/57 "2024-03-27T18:48:33Z")
My motivation for the backslash is I think trying to attract students whose only programming experience is MATLAB to Fortran would be easier if there was a way to make Fortran a little closer to what they already know from MATLAB. The interactive capabilities of LFORTRAN will go a long way in that direction. Enabling syntax that approximates MATLAB would also help. ---
**Author:** ![jkd2022](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jkd2022/32/1949_2.png) [@jkd2022](https://fortran-lang.discourse.group/u/jkd2022)\ **Post date:** [March 27, 2024, 6:51pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/58 "2024-03-27T18:51:36Z")
I feel that `inverse(...)` would be less ambiguous, so that you could form, for example: ```auto C=matmul(conjg(transpose(A)),inverse(B)) ``` ---
**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\ **Post date:** [March 27, 2024, 8:20pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/59 "2024-03-27T20:20:39Z")
> [@FedericoPerini](#): > > for those relying on `pure` BLAS/LAPACK procedures 2 out of 3 interfaces can be pure: > > - subroutine interface > - `.operator.` interface Here is how I see it… As mentioned above, one advantage of functions is that they can appear in expressions. Non-trivial expressions involving arrays can be tricky to optimize by the compilers (particularly to avoid creating too many temporaries…), and my understanding is that the `pure` or `simple` properties can help the compilers in these cases. You say that operators can do in this case: that’s correct, but style preferences set apart (I prefer functions to custom operators), not all functions will have an operator version I guess (returning the k first eigenvectors for instance require an additional argument)… maybe “most of them” is enough. Having the return flags only in the “expert” subroutine versions would look fine to me. How does Eigen manages errors? ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [March 27, 2024, 8:22pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/60 "2024-03-27T20:22:39Z")
Although the operations `A.inv.B` and `matmul(inverse(A),B)` are the same in abstract mathematics, they are very different in numerical analysis. Both the operation counts and the numerical properties differ. ---
**Author:** ![jkd2022](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jkd2022/32/1949_2.png) [@jkd2022](https://fortran-lang.discourse.group/u/jkd2022)\ **Post date:** [March 27, 2024, 8:51pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/61 "2024-03-27T20:51:13Z")
Maybe, but it’s a bit ambiguous. If I saw `A.inv.B` I would expect that A is inverted and multiplied with B rather than being the solution to _AX_ = _B_. Perhaps .div. would be better. Also, are there any provisions for declaring the matrices as symmetric or Hermitian? Finally, I think if the ‘info’ from LAPACK is not going to be returned, then if the matrix _A_ is singular the entirety of _B_ should be set to NaN. This way it fails safely. (How do you include LaTeX formulae on these forums?) ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [March 27, 2024, 10:00pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/62 "2024-03-27T22:00:22Z")
> [@jkd2022](#): > > (How do you include LaTeX formulae on these forums?) By enclosing them in a pair of dollar signs ($), looking at what others have done. The syntax is described at [LaTeX - Mathematical Python](https://patrickwalls.github.io/mathematicalpython/jupyter/latex/). [Previous page](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531.md?page=2) [Next page](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531.md?page=4)