search

ocean-transport / scale-aware-air-sea

Repo for collaborative project on scale-aware air-sea fluxes
2 stars 0 forks source link

Bulk formula package #11

Closed rabernat closed 2 years ago

rabernat commented 2 years ago

Today we discussed the choice between two options for computing bulk formulae:

For the latter, here are some helpful links:

My recommendation would be to make aerobulk a submodule of aerobulk-python.

Edited to add the following checklist

TomNicholas commented 2 years ago

Looking at f2py's documentation it seems there are a few ways of exposing the fortran routines in python. Do we want to: 1) Leave aerobulk alone entirely, make a new package "aerobulk-python" which keeps a copy of aerobulk's fortran routines in a subdirectory (via a git submodule), and create signature files using f2py which live in aerobulk-python. 2) Make PRs to aerobulk which add these F2PY directives directly into the aerobulk fortran routines as comments.

Approach (1) separates our packaging from the original routines most clearly, doesn't require the aerobulk developers to accept anything, and was what was done in xcape.

Approach (2) requires the smallest amount of additional code, and integrates the python interface more tightly with the aerobulk one, and gives us a reason to ask for commit access to aerobulk.

rabernat commented 2 years ago

I would say 1

jbusecke commented 2 years ago

Yeah, I agree. It seems that is the fastest way to make progress. If we have success with this approach, we can always consider adding code upstream later on?

TomNicholas commented 2 years ago

@rabernat I'm trying to compile the wrapper for the mod_blk_ncar.f90 file (I picked the NCAR one at random). First I tried compiling the wrapper for that file, but although I can create the signature file (mod_blk_ncar.pyf), I can't compile it with f2py immediately because it relies on a fortran module defined in mod_const.f90, and asks for a mod_const.mod file.

I don't understand how to get f2py's interface to build the .mod file for me, but I can make the .mod file by using gfortran directly on mod_const.f90.

Now using f2py on mod_blk_ncar.f90 asks for mod_phymbl.mod, which is the other module USEd by mod_blk_ncar. However if I try to compile mod_phymbl.mod with gfortran -I. then I get an error:

mod_phymbl.f90:1876:132:

 1876 |             WRITE(*,'(" *** ERROR (check_unit_consistency@mod_phymbl): shape of `mask` does not agree with array of field ",a," !")') TRIM(cfield)
      |                                                                                                                                    1
Error: Line truncated at (1) [-Werror=line-truncation]
mod_phymbl.f90:1876:132:

 1876 |             WRITE(*,'(" *** ERROR (check_unit_consistency@mod_phymbl): shape of `mask` does not agree with array of field ",a," !")') TRIM(cfield)
      |                                                                                                                                    1
Error: Syntax error in WRITE statement at (1)
mod_phymbl.f90:1949:132:

 1949 |          WRITE(*,'(" *** ERROR (check_unit_consistency@mod_phymbl): field `",a,"` does not seem to be in ",a," !")') TRIM(cfield), TRIM(cunit)
      |                                                                                                                                    1
Error: Line truncated at (1) [-Werror=line-truncation]
mod_phymbl.f90:1949:132:

 1949 |          WRITE(*,'(" *** ERROR (check_unit_consistency@mod_phymbl): field `",a,"` does not seem to be in ",a," !")') TRIM(cfield), TRIM(cunit)
      |                                                                                                                                    1
Error: Symbol ‘t’ at (1) has no IMPLICIT type
f951: some warnings being treated as errors

Am I doing this the right way? I feel like I'm doing something silly with respect to the linking of these files... Xcape has just .f90 and the .pyf signature files in its repo - does that mean that the compilation happens when the user installs it from conda?

rabernat commented 2 years ago

Have you tried to first compile the fortran as documented in the repo--indepndent of python? That will teach you how the code "expects" to be compiled and is probably a prerequisite for using f2py.

TomNicholas commented 2 years ago

Have you tried to first compile the fortran as documented in the repo--indepndent of python? That will teach you how the code "expects" to be compiled and is probably a prerequisite for using f2py.

Yeah I mean I've successfully used make with the makefile to get all the .o and .mod files. I'll try compiling and running the tests mentioned in the readme too.

I guess I'm just confused as to how to tell f2py to compile the code in the same way that the makefile specifies, or whether I can tell f2py just to use the object files I just compiled using make. (The f2py documentation curtly dismisses make in general so that isn't especially helpful).

rabernat commented 2 years ago

Some notes summarizing the progress we have made so far. Here are the steps to try to compile f2py on aerobulk. I assume you have checked out https://github.com/brodeau/aerobulk and are in the root directory.

Create a conda environment

mamba create -n aerobulk python=3.9 numpy gfortran
conda activate aerobulk

Build the aerobulk source

cp arch/make.macro_GnuLinux make.macro

those of us on mac needed to edit make.macro make this change

# These are needed for the C/C++ interface
- FF += -std=f2008 -lstdc++
+ FF += -std=f2008

EDIT:

Note that the line uncommented by default is actually using -std=gnu here. Both versions do however compile. Not sure about differences.

Then you should be able to run make

This will build:

I was able to run one of the binaries with bin/example_call_aerobulk.x

output ``` *********** COARE 3.0 ***************** =================================================================== ----- AeroBulk_init ----- *** Bulk parameterization to be used => "coare3p0" ==> will use the Cool-skin & Warm-ayer scheme of `coare3p0` ! *** Computational domain shape: Ni x Nj = 00002 x 00001 *** Number of time records that will be treated: 1 *** Number of iterations in bulk algos: nb_iter = 10 *** Filling the `mask` array... ==> no points need to be masked! :) *** Type of prescribed air humidity `specific humidity [kg/kg]` =================================================================== =================================================================== ----- AeroBulk_bye ----- =================================================================== --------------------------------------------------------------------- Parameter | Unstable ASL | Stable ASL | units --------------------------------------------------------------------- Wind speed at zu = 5.00000000 5.00000000 m/s SST = 22.0000000 22.0000000 deg.C Abs. temperature at zt = 20.0000000 25.0000000 deg.C Pot. temperature at zt = 20.0134144 25.0150242 deg.C Sensible heat flux: QH = -15.1545086 17.8423691 W/m**2 Latent heat flux: QL = -81.3846741 -50.8408585 W/m**2 Evaporation: Evap = -2.87061906 -1.79333246 mm/day Skin temperature: SSST = 21.7219677 21.7576351 deg.C Tau_x = 3.57834995E-02 1.73626728E-02 N/m**2 Tau_y = 0.00000000 0.00000000 N/m**2 Tau = 3.57834995E-02 1.73626728E-02 N/m**2 *********** COARE 3.6 ***************** =================================================================== ----- AeroBulk_init ----- *** Bulk parameterization to be used => "coare3p6" ==> will use the Cool-skin & Warm-ayer scheme of `coare3p6` ! *** Computational domain shape: Ni x Nj = 00002 x 00001 *** Number of time records that will be treated: 1 *** Number of iterations in bulk algos: nb_iter = 10 *** Filling the `mask` array... ==> no points need to be masked! :) *** Type of prescribed air humidity `specific humidity [kg/kg]` =================================================================== =================================================================== ----- AeroBulk_bye ----- =================================================================== --------------------------------------------------------------------- Parameter | Unstable ASL | Stable ASL | units --------------------------------------------------------------------- Wind speed at zu = 5.00000000 5.00000000 m/s SST = 22.0000000 22.0000000 deg.C Abs. temperature at zt = 20.0000000 25.0000000 deg.C Pot. temperature at zt = 20.0134144 25.0150242 deg.C Sensible heat flux: QH = -15.3865499 17.0818920 W/m**2 Latent heat flux: QL = -83.0788422 -48.4459152 W/m**2 Evaporation: Evap = -2.93033028 -1.70883954 mm/day Skin temperature: SSST = 21.7057972 21.7485600 deg.C Tau_x = 3.21817845E-02 1.51577257E-02 N/m**2 Tau_y = 0.00000000 0.00000000 N/m**2 Tau = 3.21817845E-02 1.51577257E-02 N/m**2 *********** ECMWF ***************** =================================================================== ----- AeroBulk_init ----- *** Bulk parameterization to be used => "ecmwf" ==> will use the Cool-skin & Warm-ayer scheme of `ecmwf` ! *** Computational domain shape: Ni x Nj = 00002 x 00001 *** Number of time records that will be treated: 1 *** Number of iterations in bulk algos: nb_iter = 10 *** Filling the `mask` array... ==> no points need to be masked! :) *** Type of prescribed air humidity `specific humidity [kg/kg]` =================================================================== =================================================================== ----- AeroBulk_bye ----- =================================================================== --------------------------------------------------------------------- Parameter | Unstable ASL | Stable ASL | units --------------------------------------------------------------------- Wind speed at zu = 5.00000000 5.00000000 m/s SST = 22.0000000 22.0000000 deg.C Abs. temperature at zt = 20.0000000 25.0000000 deg.C Pot. temperature at zt = 20.0134144 25.0150242 deg.C Sensible heat flux: QH = -14.3822346 17.6531811 W/m**2 Latent heat flux: QL = -80.2958984 -52.4623947 W/m**2 Evaporation: Evap = -2.83224440 -1.85053933 mm/day Skin temperature: SSST = 21.7325401 21.7630310 deg.C Tau_x = 3.84389125E-02 1.93256922E-02 N/m**2 Tau_y = 0.00000000 0.00000000 N/m**2 Tau = 3.84389125E-02 1.93256922E-02 N/m**2 *********** NCAR ***************** =================================================================== ----- AeroBulk_init ----- *** Bulk parameterization to be used => "ncar" *** Cool-skin & Warm-layer schemes will NOT be used! *** Computational domain shape: Ni x Nj = 00002 x 00001 *** Number of time records that will be treated: 1 *** Number of iterations in bulk algos: nb_iter = 10 *** Filling the `mask` array... ==> no points need to be masked! :) *** Type of prescribed air humidity `specific humidity [kg/kg]` =================================================================== =================================================================== ----- AeroBulk_bye ----- =================================================================== --------------------------------------------------------------------- Parameter | Unstable ASL | Stable ASL | units --------------------------------------------------------------------- Wind speed at zu = 5.00000000 5.00000000 m/s SST = 22.0000000 22.0000000 deg.C Abs. temperature at zt = 20.0000000 25.0000000 deg.C Pot. temperature at zt = 20.0134144 25.0150242 deg.C Sensible heat flux: QH = -16.6969528 10.7261686 W/m**2 Latent heat flux: QL = -88.4781876 -71.9012222 W/m**2 Evaporation: Evap = -3.12166286 -2.53679895 mm/day Tau_x = 3.58519591E-02 2.77329944E-02 N/m**2 Tau_y = 0.00000000 0.00000000 N/m**2 Tau = 3.58519591E-02 2.77329944E-02 N/m**2 *********** ANDREAS ***************** =================================================================== ----- AeroBulk_init ----- *** Bulk parameterization to be used => "andreas" *** Cool-skin & Warm-layer schemes will NOT be used! *** Computational domain shape: Ni x Nj = 00002 x 00001 *** Number of time records that will be treated: 1 *** Number of iterations in bulk algos: nb_iter = 10 *** Filling the `mask` array... ==> no points need to be masked! :) *** Type of prescribed air humidity `specific humidity [kg/kg]` =================================================================== =================================================================== ----- AeroBulk_bye ----- =================================================================== --------------------------------------------------------------------- Parameter | Unstable ASL | Stable ASL | units --------------------------------------------------------------------- Wind speed at zu = 5.00000000 5.00000000 m/s SST = 22.0000000 22.0000000 deg.C Abs. temperature at zt = 20.0000000 25.0000000 deg.C Pot. temperature at zt = 20.0134144 25.0150242 deg.C Sensible heat flux: QH = -14.4129944 15.1970539 W/m**2 Latent heat flux: QL = -74.4637756 -51.7018318 W/m**2 Evaporation: Evap = -2.62721038 -1.82412970 mm/day Tau_x = 3.02770734E-02 1.79436151E-02 N/m**2 Tau_y = 0.00000000 0.00000000 N/m**2 Tau = 3.02770734E-02 1.79436151E-02 N/m**2 ```

Do the f2py magic

First we build the signature file

cd src
python -m numpy.f2py mod_aerobulk.f90 -m mod_aerobulk -h mod_aerobulk.pyf
mod_aerobulk.pyf ``` ! -*- f90 -*- ! Note: the context of this file is case sensitive. python module mod_aerobulk ! in interface ! in :mod_aerobulk module mod_aerobulk ! in :mod_aerobulk:mod_aerobulk.f90 use mod_const use mod_phymbl, only: type_of_humidity,check_unit_consistency use mod_aerobulk_compute subroutine aerobulk_init(nt,calgo,psst,pta,pha,pu,pv,pslp,l_use_skin,prsw,prlw) ! in :mod_aerobulk:mod_aerobulk.f90:mod_aerobulk integer intent(in) :: nt character*(*) intent(in) :: calgo real(kind=wp) dimension(:,:),intent(in) :: psst real(kind=wp) dimension(:,:),intent(in) :: pta real(kind=wp) dimension(:,:),intent(in) :: pha real(kind=wp) dimension(:,:),intent(in) :: pu real(kind=wp) dimension(:,:),intent(in) :: pv real(kind=wp) dimension(:,:),intent(in) :: pslp logical, optional,intent(in) :: l_use_skin real(kind=wp), optional,dimension(:,:),intent(in) :: prsw real(kind=wp), optional,dimension(:,:),intent(in) :: prlw end subroutine aerobulk_init subroutine aerobulk_bye ! in :mod_aerobulk:mod_aerobulk.f90:mod_aerobulk end subroutine aerobulk_bye subroutine aerobulk_model(jt,nt,calgo,zt,zu,sst,t_zt,hum_zt,u_zu,v_zu,slp,ql,qh,tau_x,tau_y,evap,niter,l_use_skin,rad_sw,rad_lw,t_s) ! in :mod_aerobulk:mod_aerobulk.f90:mod_aerobulk integer intent(in) :: jt integer intent(in) :: nt character*(*) intent(in) :: calgo real(kind=wp) intent(in) :: zt real(kind=wp) intent(in) :: zu real(kind=wp) dimension(:,:),intent(in) :: sst real(kind=wp) dimension(:,:),intent(in) :: t_zt real(kind=wp) dimension(:,:),intent(in) :: hum_zt real(kind=wp) dimension(:,:),intent(in) :: u_zu real(kind=wp) dimension(:,:),intent(in) :: v_zu real(kind=wp) dimension(:,:),intent(in) :: slp real(kind=wp) dimension(:,:),intent(out) :: ql real(kind=wp) dimension(:,:),intent(out) :: qh real(kind=wp) dimension(:,:),intent(out) :: tau_x real(kind=wp) dimension(:,:),intent(out) :: tau_y real(kind=wp) dimension(:,:),intent(out) :: evap integer, optional,intent(in) :: niter logical, optional,intent(in) :: l_use_skin real(kind=wp), optional,dimension(:,:),intent(in) :: rad_sw real(kind=wp), optional,dimension(:,:),intent(in) :: rad_lw real(kind=wp), optional,dimension(:,:),intent(out) :: t_s end subroutine aerobulk_model end module mod_aerobulk end interface end python module mod_aerobulk ! This file was auto-generated with f2py (version:1.22.3). ! See: ! https://web.archive.org/web/20140822061353/http://cens.ioc.ee/projects/f2py2e ```

At this point, we realized that, since the aerobulk code uses KIND specifiers to define fortran types, we need to follow these instructions - https://numpy.org/doc/stable/f2py/advanced.html#dealing-with-kind-specifiers

We created a file src/.f2py_f2cmap with the contents

{'real': {'wp': 'double', 'dp': 'double', 'sp': 'float'}}

Then the magic command to build

python -m numpy.f2py --verbose -c --f90flags="-fdefault-real-8 -ffree-line-length-200 --std=gnu" \
   -I../mod mod_aerobulk.pyf mod_*.f90

Unfortunately, this is still failing with the error

/var/folders/n8/63q49ms55wxcj_gfbtykwp5r0000gn/T/tmp1gwaq670/src.macosx-10.9-x86_64-3.9/mod_aerobulk-f2pywrappers2.f90:30:16:

   30 |       real(kind=wp) psst(f2py_psst_d0,f2py_psst_d1)
      |                1
Error: Parameter ‘wp’ at (1) has not been declared or is a variable, which does not reduce to a constant expression

That's where we are at.

rabernat commented 2 years ago

Ok, quick update, replacing kind=wp with kind=8 in mod_aerobulk.pyf allows the f2py to succeed. :tada:

rabernat commented 2 years ago

Very productive session yesterday. I updated the top level issue with a checklist of next steps.

@jbusecke - since you are our packaging expert, I think it makes sense for you to create the repo and set up the package skeleton. Shall we do this in the xgcm org?

In the meanwhile, maybe @TomNicholas can keep tweaking f2py to address items 2 and 3?

Does this sounds like a good plan?

TomNicholas commented 2 years ago

Ok, quick update, replacing kind=wp with kind=8 in mod_aerobulk.pyf allows the f2py to succeed.

I'm having problems replicating this step on my machine. After doing the substitution, I still get compilation errors of the form

x86_64-conda-linux-gnu-gfortran:f90: mod_blk_ncar-f2pywrappers2.f90
mod_blk_ncar-f2pywrappers2.f90:16:16:

   16 |       real(kind=wp) zt
      |                1
Error: Parameter 'wp' at (1) has not been declared or is a variable, which does not reduce to a constant expression

These are similar to what we had before changing kind=wp to kind=8 for mod_aerobulk, except now it's happening with a file that I shouldn't need a signature for mod_blk_ncar.

rabernat commented 2 years ago

After doing the substitution, I still get compilation errors

Did you remember to add the .f2py_f2cmap file?

TomNicholas commented 2 years ago

Did you remember to add the .f2py_f2cmap file?

Yes, and the compilation prints 

Reading f2cmap from '.f2py_f2cmap' ...
    Mapping "real(kind=wp)" to "double"
    Mapping "real(kind=dp)" to "double"
    Mapping "real(kind=sp)" to "float"
Successfully applied user defined f2cmap changes
rabernat commented 2 years ago

So the question is, what process is generating the mod_blk_ncar-f2pywrappers2.f90 file? Presumably a manual search / replace within that file (and all others like it) would resolve the issue.

jbusecke commented 2 years ago

Closing now that we have all the intstructions in https://github.com/xgcm/aerobulk-python/pull/19