From a551f750e3f608916f852a3bcce90efce162f62d Mon Sep 17 00:00:00 2001 From: Steve Jester Date: Fri, 20 Sep 2024 14:01:09 -0400 Subject: [PATCH 1/6] created and working on wrf_ghg.f90 for use with getvar --- fortran/wrf_ghg.f90 | 44 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 fortran/wrf_ghg.f90 diff --git a/fortran/wrf_ghg.f90 b/fortran/wrf_ghg.f90 new file mode 100644 index 0000000..5131828 --- /dev/null +++ b/fortran/wrf_ghg.f90 @@ -0,0 +1,44 @@ +SUBROUTINE DCOMPUTEXGHG(sfc_p, pres, nx, ny, nz, ant, bio, bck, xghg) + IMPLICIT NONE + + !f2py threadsafe + !f2py intent(in,out) :: xghg + + INTEGER, INTENT(IN) :: nx, ny, nz + REAL(KIND=8), DIMENSION(nx, ny), INTENT(IN) :: sfc_p + REAL(KIND=8), DIMENSION(nx, ny, nz), INTENT(IN) :: pres + REAL(KIND=8), DIMENSION(nx, ny, nz), INTENT(IN) :: ant + REAL(KIND=8), DIMENSION(nx, ny, nz), INTENT(IN) :: bio + REAL(KIND=8), DIMENSION(nx, ny, nz), INTENT(IN) :: bck + REAL(KIND=8), DIMENSION(nx, ny), INTENT(IN,OUT) :: xghg + + INTEGER :: i, j, k + REAL(KIND=8), DIMENSION(nx, ny, nz) :: ghg + REAL(KIND=8), DIMENSION(nx, ny, nz) :: pres_bound + REAL(KIND=8), DIMENSION(nx, ny, nz) :: p_layer_diff + REAL(KIND=8), DIMENSION(nx, ny) :: p_diff + + + + !$OMP PARALLEL + !$OMP DO COLLAPSE(3) SCHEDULE(runtime) + DO k = 1,nz + DO j = 1,ny + DO i = 1,nx + ghg(i , j, k) = ant(i, j, k) + bio(i, j, k) - bck(i, j, k) + IF (k .eq. 1) THEN + pres_bound(i, j, k) = sfc_p(i, j) + ELSE + pres_bound(i, j, k) = pres_bound(i, j, k-1) + (2*(pres(i, j, k-1)-pres_bound(i, j, k-1))) + END IF + END DO + END DO + END DO + !$OMP END DO + + !$OMP DO COLLAPSE(2) SCHEDULE(runtime) + DO j = i,ny + DO i = 1,nx + END DO + END DO + From 35513d69728fd86b1e00d7a9d5f680d83f15162f Mon Sep 17 00:00:00 2001 From: Steve Jester Date: Sat, 21 Sep 2024 12:18:15 -0400 Subject: [PATCH 2/6] completed ver 1 of the get_ghg and related functions --- fortran/wrf_ghg.f90 | 13 +++++- setup.py | 3 +- src/wrf/extension.py | 24 +++++++++- src/wrf/g_ghg.py | 102 +++++++++++++++++++++++++++++++++++++++++++ src/wrf/routines.py | 9 +++- 5 files changed, 146 insertions(+), 5 deletions(-) create mode 100644 src/wrf/g_ghg.py diff --git a/fortran/wrf_ghg.f90 b/fortran/wrf_ghg.f90 index 5131828..39ef9cb 100644 --- a/fortran/wrf_ghg.f90 +++ b/fortran/wrf_ghg.f90 @@ -17,6 +17,7 @@ SUBROUTINE DCOMPUTEXGHG(sfc_p, pres, nx, ny, nz, ant, bio, bck, xghg) REAL(KIND=8), DIMENSION(nx, ny, nz) :: pres_bound REAL(KIND=8), DIMENSION(nx, ny, nz) :: p_layer_diff REAL(KIND=8), DIMENSION(nx, ny) :: p_diff + REAL(KIND=8), DIMENSION(nx, ny) :: weighted_ghg @@ -28,17 +29,27 @@ SUBROUTINE DCOMPUTEXGHG(sfc_p, pres, nx, ny, nz, ant, bio, bck, xghg) ghg(i , j, k) = ant(i, j, k) + bio(i, j, k) - bck(i, j, k) IF (k .eq. 1) THEN pres_bound(i, j, k) = sfc_p(i, j) + p_layer_diff(i, j, k) = pres(i, j, k) - sfc_p(i, j) ELSE pres_bound(i, j, k) = pres_bound(i, j, k-1) + (2*(pres(i, j, k-1)-pres_bound(i, j, k-1))) + p_layer_diff(i, j, k) = pres_bound(i, j, k-1) - pres_bound(i, j, k) END IF END DO END DO END DO !$OMP END DO - + + weighted_ghg = SUM(ghg * p_layer_diff, 3) + !$OMP DO COLLAPSE(2) SCHEDULE(runtime) DO j = i,ny DO i = 1,nx + p_diff(i, j) = pres_bound(i, j, 1) - pres_bound(i, j, nz) + xghg(i, j) = weighted_ghg(i, j) / p_diff(i, j) END DO END DO + !$OMP END DO + !$OMP END PARALLEL + RETURN +END SUBROUTINE DCOMPUTEXGHG \ No newline at end of file diff --git a/setup.py b/setup.py index 0bbf359..deddabd 100755 --- a/setup.py +++ b/setup.py @@ -50,7 +50,8 @@ "fortran/wrf_pw.f90", "fortran/wrf_vinterp.f90", "fortran/wrf_wind.f90", - "fortran/omp.f90"] + "fortran/omp.f90" + "fortran/wrf_ghg.f90"] ) #Note: __version__ will be set in the version.py script loaded below diff --git a/src/wrf/extension.py b/src/wrf/extension.py index d28d758..70ddc0c 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -12,7 +12,7 @@ dlltoij, dijtoll, deqthecalc, omgcalc, virtual_temp, wetbulbcalc, dcomputepw, wrf_monotonic, wrf_vintrp, dcomputewspd, - dcomputewdir, dinterp3dz_2dlev, + dcomputewdir, dinterp3dz_2dlev, dcomputexghg, fomp_set_num_threads, fomp_get_num_threads, fomp_get_max_threads, fomp_get_thread_num, fomp_get_num_procs, fomp_in_parallel, @@ -439,6 +439,28 @@ def _uvmet(u, v, lat, lon, cen_long, cone, isstag=0, has_missing=False, return result +@check_args(0, 2, (2, 3, 3, 3, 3)) +@left_iteration(2, 2, ref_var_idx=0) +@cast_type(arg_idxs=(0, 1, 2, 3, 4)) +@extract_and_transpose() +def _xghg(sfc_p, pres, ant, bio, bck, outview=None): + """ + Wrapper for dcomputexghg. + + Located in wrf_ghg.f90 + """ + if outview is None: + outview = np.empty_like(sfc_p) + + result = dcomputexghg(sfc_p, + pres, + ant, + bio, + bck, + outview) + + return result + @check_args(0, 3, (3, 3, 3, 3)) @left_iteration(3, 3, ref_var_idx=0) diff --git a/src/wrf/g_ghg.py b/src/wrf/g_ghg.py new file mode 100644 index 0000000..d93df1b --- /dev/null +++ b/src/wrf/g_ghg.py @@ -0,0 +1,102 @@ +from __future__ import (absolute_import, division, print_function) +import numpy as np +from .extension import _xghg +from .metadecorators import copy_and_set_metadata +from .util import extract_vars + +@copy_and_set_metadata(copy_varname='CH4_BCK', name='column_averaged_ghg', + description='column averaged ghg', + units='ppmv') +def get_xghg(wrfin, ghg="xch4", timeidx=0, method="cat", squeeze=True, + cache=None, meta=True, _key=None): + """Return column averaged greenhouse gas mixing ratio. + + This functions extracts the necessary variables from the NetCDF file + object in order to perform the calculation. + + Args: + + wrfin (:class:`netCDF4.Dataset`, :class:`Nio.NioFile`, or an \ + iterable): WRF-ARW NetCDF + data as a :class:`netCDF4.Dataset`, :class:`Nio.NioFile` + or an iterable sequence of the aforementioned types. + + ghg (:obj:`str`, optional): The desired greenhouse gas. Must be + one of 'xch4', 'xco2', or 'xco'. The default is 'xch4'. A + non-supported choice raises a NotImplementedError. + + timeidx (:obj:`int` or :data:`wrf.ALL_TIMES`, optional): The + desired time index. This value can be a positive integer, + negative integer, or + :data:`wrf.ALL_TIMES` (an alias for None) to return + all times in the file or sequence. The default is 0. + + method (:obj:`str`, optional): The aggregation method to use for + sequences. Must be either 'cat' or 'join'. + 'cat' combines the data along the Time dimension. + 'join' creates a new dimension for the file index. + The default is 'cat'. + + squeeze (:obj:`bool`, optional): Set to False to prevent dimensions + with a size of 1 from being automatically removed from the + shape of the output. Default is True. + + cache (:obj:`dict`, optional): A dictionary of (varname, ndarray) + that can be used to supply pre-extracted NetCDF variables to + the computational routines. It is primarily used for internal + purposes, but can also be used to improve performance by + eliminating the need to repeatedly extract the same variables + used in multiple diagnostics calculations, particularly when + using large sequences of files. + Default is None. + + meta (:obj:`bool`, optional): Set to False to disable metadata and + return :class:`numpy.ndarray` instead of + :class:`xarray.DataArray`. Default is True. + + _key (:obj:`int`, optional): A caching key. This is used for + internal purposes only. Default is None. + + Returns: + :class:`xarray.DataArray` or :class:`numpy.ndarray`: Omega. + If xarray is + enabled and the *meta* parameter is True, then the result will be a + :class:`xarray.DataArray` object. Otherwise, the result will be a + :class:`numpy.ndarray` object with no metadata. + """ + + match ghg: + case 'xch4' | 'ch4' | 'methane': + varnames = ('P', 'PSFC', 'CH4_ANT', 'CH4_BIO','CH4_BCK') + ncvars = extract_vars(wrfin, timeidx, varnames, method, squeeze, cache, + meta=False, _key=_key) + pres = ncvars['P'] + sfc_p = ncvars['PSFC'] + ant = ncvars['CH4_ANT'] + bio = ncvars['CH4_BIO'] + bck = ncvars['CH4_BCK'] + case 'xco2' | 'co2' | 'carbon dioxide': + varnames = ('P', 'PSFC', 'CO2_ANT', 'CO2_BIO','CO2_BCK') + ncvars = extract_vars(wrfin, timeidx, varnames, method, squeeze, cache, + meta=False, _key=_key) + pres = ncvars['P'] + sfc_p = ncvars['PSFC'] + ant = ncvars['CO2_ANT'] + bio = ncvars['CO2_BIO'] + bck = ncvars['CO2_BCK'] + case 'xco' | 'co' | 'carbon monoxide': + varnames = ('P', 'PSFC', 'CO_ANT', 'CO_BCK') + ncvars = extract_vars(wrfin, timeidx, varnames, method, squeeze, cache, + meta=False, _key=_key) + pres = ncvars['P'] + sfc_p = ncvars['PSFC'] + ant = ncvars['CO_ANT'] + bio = np.zeros_like(ant) + bck = ncvars['CO_BCK'] + + case _: + raise NotImplementedError('Chemistry variable not supported at this time.') + + xghg = _xghg(sfc_p, pres, ant, bio, bck) + + return xghg \ No newline at end of file diff --git a/src/wrf/routines.py b/src/wrf/routines.py index d40b72d..1230f6e 100644 --- a/src/wrf/routines.py +++ b/src/wrf/routines.py @@ -31,6 +31,7 @@ from .g_times import get_times, get_xtimes from .g_cloudfrac import (get_cloudfrac, get_low_cloudfrac, get_mid_cloudfrac, get_high_cloudfrac) +from .g_ghg import get_xghg # func is the function to call. kargs are required arguments that should @@ -97,7 +98,8 @@ "wdir10": get_destag_wdir10, "low_cloudfrac": get_low_cloudfrac, "mid_cloudfrac": get_mid_cloudfrac, - "high_cloudfrac": get_high_cloudfrac + "high_cloudfrac": get_high_cloudfrac, + "column_averaged_ghg": get_xghg } _VALID_KARGS = {"cape2d": ["missing"], @@ -166,6 +168,7 @@ "mid_thresh", "high_thresh"], "high_cloudfrac": ["vert_type", "low_thresh", "mid_thresh", "high_thresh"], + "column_averaged_ghg": [], "default": [] } @@ -197,7 +200,9 @@ "wspd_uvmet10": "uvmet10_wspd", "wdir_uvmet10": "uvmet10_wdir", "mcape": "cape2d_only", - "mcin": "cin2d_only" + "mcin": "cin2d_only", + "xghg": "column_averaged_ghg", + "ghg": "column_averaged_ghg" } From 1a86a356e481a6c638cb4b50c330d3317d85c48a Mon Sep 17 00:00:00 2001 From: Steve Jester Date: Sat, 21 Sep 2024 12:23:10 -0400 Subject: [PATCH 3/6] added ghg valid kwarg --- src/wrf/g_ghg.py | 4 ++-- src/wrf/routines.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/wrf/g_ghg.py b/src/wrf/g_ghg.py index d93df1b..ba78465 100644 --- a/src/wrf/g_ghg.py +++ b/src/wrf/g_ghg.py @@ -7,8 +7,8 @@ @copy_and_set_metadata(copy_varname='CH4_BCK', name='column_averaged_ghg', description='column averaged ghg', units='ppmv') -def get_xghg(wrfin, ghg="xch4", timeidx=0, method="cat", squeeze=True, - cache=None, meta=True, _key=None): +def get_xghg(wrfin, timeidx=0, method="cat", squeeze=True, + cache=None, meta=True, _key=None, ghg="xch4"): """Return column averaged greenhouse gas mixing ratio. This functions extracts the necessary variables from the NetCDF file diff --git a/src/wrf/routines.py b/src/wrf/routines.py index 1230f6e..60b4bb4 100644 --- a/src/wrf/routines.py +++ b/src/wrf/routines.py @@ -168,7 +168,7 @@ "mid_thresh", "high_thresh"], "high_cloudfrac": ["vert_type", "low_thresh", "mid_thresh", "high_thresh"], - "column_averaged_ghg": [], + "column_averaged_ghg": ["ghg"], "default": [] } From 8759e97de44413b8e87687ce0db752aee40319e9 Mon Sep 17 00:00:00 2001 From: Steve Jester Date: Sat, 21 Sep 2024 16:41:08 -0400 Subject: [PATCH 4/6] forgot a comma in setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index deddabd..2ff5bb0 100755 --- a/setup.py +++ b/setup.py @@ -50,7 +50,7 @@ "fortran/wrf_pw.f90", "fortran/wrf_vinterp.f90", "fortran/wrf_wind.f90", - "fortran/omp.f90" + "fortran/omp.f90", "fortran/wrf_ghg.f90"] ) From 5c7834ca167737830ea17947584daa893bbf566b Mon Sep 17 00:00:00 2001 From: Steve Jester Date: Sat, 21 Sep 2024 16:58:55 -0400 Subject: [PATCH 5/6] changed the intent of xghg in wrf_ghg.f90 --- fortran/wrf_ghg.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fortran/wrf_ghg.f90 b/fortran/wrf_ghg.f90 index 39ef9cb..678ade6 100644 --- a/fortran/wrf_ghg.f90 +++ b/fortran/wrf_ghg.f90 @@ -10,7 +10,7 @@ SUBROUTINE DCOMPUTEXGHG(sfc_p, pres, nx, ny, nz, ant, bio, bck, xghg) REAL(KIND=8), DIMENSION(nx, ny, nz), INTENT(IN) :: ant REAL(KIND=8), DIMENSION(nx, ny, nz), INTENT(IN) :: bio REAL(KIND=8), DIMENSION(nx, ny, nz), INTENT(IN) :: bck - REAL(KIND=8), DIMENSION(nx, ny), INTENT(IN,OUT) :: xghg + REAL(KIND=8), DIMENSION(nx, ny), INTENT(OUT) :: xghg INTEGER :: i, j, k REAL(KIND=8), DIMENSION(nx, ny, nz) :: ghg From fbe0ba2fc109300174f2923a6ab188ec5b14b99c Mon Sep 17 00:00:00 2001 From: Steve Jester Date: Thu, 3 Oct 2024 12:46:30 -0400 Subject: [PATCH 6/6] added test, initial test failed, recompiling --- .gitignore | 1 + fortran/build_help/sizes | Bin 0 -> 16432 bytes fortran/omp.f90 | 332 ++++++++++++++++--- fortran/ompgen.F90 | 2 +- wrf_ghg_test.ipynb | 674 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 965 insertions(+), 44 deletions(-) create mode 100755 fortran/build_help/sizes create mode 100644 wrf_ghg_test.ipynb diff --git a/.gitignore b/.gitignore index a40dae0..79920cf 100644 --- a/.gitignore +++ b/.gitignore @@ -56,3 +56,4 @@ build .settings src/wrf_python.egg-info +wrfout_d02_2023-07-21_18:00:00 diff --git a/fortran/build_help/sizes b/fortran/build_help/sizes new file mode 100755 index 0000000000000000000000000000000000000000..8f2ba95d8cd49878a2c8e7208a991260a6233a31 GIT binary patch literal 16432 zcmeGjYiu0FdCnITQ^MJSpqfy^HceV5s?T;38>c9+=U3LP^QdD(N&}no-TLm}9?sof zV)GE8#;8YdYqf%^54C~vBciB6RcUEd#KTbtk4Er=mQvJ46;w7!R8s;aJhJ`1nfY$_ z)^~N=AN0?Tb^Fctp80NOZpJ&`9o>CBRX!i1;%8Sg#2&2^6jd@FZI=}QDr;vKz;gw= zoSg^wB2E*0RRG|M60a(SHInWJh5u#qDa%{2S5M;6+s8@|FiM!$L zoXHqL;urlXs0>Bn&mugA0#FD#^)NO&8#$Zp+*AtZWWV2+kn|8my|~niOFe?Sr9453 zC&mOnMg|#YH|&z~1W^{tS3wWk`JWaFd%yHoDK74DAj8X= zl8K?V4Nb{tV=|G>k2a3Bw>7qH2<9@ub-dqbA7wXh9pp`-I4T@t49D9N_~7vP6nrj# zkF%gI>QMhs9)|?R_k8$JUr`2iQC`sQgqm&dW#f2?6Dm9|2f9**cN%CHm7(`kpre1K zjQkfX(C+|xIs35+^am@@U#LKTyaN5tKwk-;c}4FX*rG>`tPxA(Od~t6r8AjH8w25? zq``DOmdd2{oEgrVy3V+0dBNClCTnKHX&oQAVI!+2(xwpuxG$S9!BSBmX9|)Y1vd_Q zQ85w^XTuR_%6wJB)H5SyB9qQJ6v?2I1cKpg45akPXjmUkq{GR?E!ZCzy%E`;R5+1l zfK$tG47Va7g|#-$%A`gRY2io&G>h1td^nm4n{gyGBYVV*XN_>Qh)1y2#y0nC>mTS3 zZPnY^=DyyJPJLZ)9k^+VXHrHJ1bna2luD!%kxV)oUf*P-_vV`RWry{VF*6R4AIc|^ zQJiw{wfSIP_~Ge;C;e4%+yXAnI|V-Rv)_TK>5GGxhPF@!b=W1&1q!3>Y&RdpUavGNT0x~w6*XLEXw5)My^2q<{veNvqKX>QsBcJVQ1L0`as8pvtcn`S<2vRQ zF9Kc!ya;#^@FL(vz>C0t8v$+nRE>6^`iUmSw27z8>X`#zK^s3(Gc7DxtgDq|H!u1?%ItiX; zt^>t0+T^AfsA!YbyAf=B-wZ%UUk^eZI$C9&mF^k!TwD(pv3YGgWC9}`+&)`!kgw2?ZBp4KxX=&2{U8xI9h!I zaNiBn)2Mawz#xcecb=*PZ_u>K0P0P=|FLB~JHgIb76#(ki8;Xc^V%7xVKB5wPSgs| z&ioO?BzhRdq4NXW%yH!2YI702jic48-G0#K0)F<`mJZ|h1%Ccy5W z*Y+0P3k~)bUJ1d9vRZ9mY*3qQM7X1G&6_QM(w_d>OuPMkJZ)LwBHT3~&J zxB3e2b%70Q+3VW)v`^c3GXEAXpkKcsv@3K&=*EyfJvAI8`FCJl)TQ1|Vs?m+#w{M`X%a#44n;ofR3(0uoK zTA+P=@s@y^477&=&7nX;M?mQa)PsCSpoYIsdjr%T2Ro(T92)hC7XdE#zdq%^?o>ov-TV0*AtS=Zds+T6Boz0%NcL=`P;3VP%Ejcb?? z=66pU;R0{WP-2mY(y}4g5^Qd4StD!YqAbK=vWUG+)2S$f9^o=ei+KQ~@&~B+VD`{F zBl>Axxp-b2&Ph70BP72R8tP-0F!wlo9_R|IXYDFU3W)3AMG%J(Iq_*;9hY>mw^zjPze**NSxTq?^2CaiXuON?kB6Vf@`wGy@-Ou+! z$5-0TE7n#3eQ8yFaeNYEt)#0eQ3Sja=onA<3V7lFgM!*-P8T%i^CqCzImSQn_EwNT z3j!#=T_s5YzgI#2F>Zf;`+r+O{xwbyuqWlhPHsY+4!;(zdCekOGiT<9hl3GjXBp^b zN{?Wk0bbFMX7pGxGZaqhQ8Sax>EZk+i@^UslZI(T;c|Uh3B1&wfJ^?_@R)9-&FmN( z&W2Nl9?hpxW1!-|AQd5j`2eSy}qzHJ?$K7CtH&yMZ^eIV4)*NuWDuA-#`o<_m+uCHwkZRza<9^C7X=#gBW z_XY6MWD9q$g6hm4u=6R3c@rfG7pf!g;NNC=^lpQ6h-f0M=W_o2Td3r5G%pe&%76`Z|BHV64FC5R5ZtdR_(f~<`sO%~*V2|_g(%K&N`qwvIO z0YVwxfuIqW6CoZ2BZLW6F(ZTq)xkqJm56}*jER<^aTtgoM3x0%PNrb0m98uQD~9`8 zSZd)nRup%>8U3zGHkYM3!q346_l~4rDbEgNvC!{M{pC--p;4?H(UH-}5luMuXWd!NmMe;6l?0%_F-yc(? zq=dU9O2FNp9|HsXPx`wQAxW@Hpv9-#|0BS0>K{~uD8U9~Its;&NVqS@vk1;1EVNIi zd|0@P0k" +# 1 "" +# 1 "ompgen.F90" MODULE omp_constants + + USE omp_lib + INTEGER, PARAMETER :: fomp_sched_kind = 4 INTEGER, PARAMETER :: fomp_lock_kind = 4 INTEGER, PARAMETER :: fomp_nest_lock_kind = 8 - INTEGER(KIND=fomp_sched_kind), PARAMETER :: fomp_sched_static = 1 - INTEGER(KIND=fomp_sched_kind), PARAMETER :: fomp_sched_dynamic = 2 - INTEGER(KIND=fomp_sched_kind), PARAMETER :: fomp_sched_guided = 3 - INTEGER(KIND=fomp_sched_kind), PARAMETER :: fomp_sched_auto = 4 + INTEGER(KIND=4), PARAMETER :: fomp_sched_static = 1 + INTEGER(KIND=4), PARAMETER :: fomp_sched_dynamic = 2 + INTEGER(KIND=4), PARAMETER :: fomp_sched_guided = 3 + INTEGER(KIND=4), PARAMETER :: fomp_sched_auto = 4 + +# 22 "ompgen.F90" + END MODULE omp_constants @@ -18,140 +27,220 @@ FUNCTION fomp_enabled() LOGICAL :: fomp_enabled - fomp_enabled = .FALSE. + + fomp_enabled = .TRUE. + + + END FUNCTION fomp_enabled SUBROUTINE fomp_set_num_threads(num_threads) + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER, INTENT(IN) :: num_threads - IF (.FALSE.) PRINT *, num_threads + + + CALL omp_set_num_threads(num_threads) + + + + END SUBROUTINE fomp_set_num_threads FUNCTION fomp_get_num_threads() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER :: fomp_get_num_threads - fomp_get_num_threads = -1 + + fomp_get_num_threads = omp_get_num_threads() + + + END FUNCTION fomp_get_num_threads FUNCTION fomp_get_max_threads() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER :: fomp_get_max_threads - fomp_get_max_threads = -1 + + fomp_get_max_threads = omp_get_max_threads() + + + END FUNCTION fomp_get_max_threads FUNCTION fomp_get_thread_num() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER :: fomp_get_thread_num - fomp_get_thread_num = -1 + + fomp_get_thread_num = omp_get_thread_num() + + + END FUNCTION fomp_get_thread_num FUNCTION fomp_get_num_procs() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER :: fomp_get_num_procs - fomp_get_num_procs = -1 + + fomp_get_num_procs = omp_get_num_procs() + + + END FUNCTION fomp_get_num_procs FUNCTION fomp_in_parallel() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe LOGICAL :: fomp_in_parallel - fomp_in_parallel = .FALSE. + + fomp_in_parallel = omp_in_parallel() + + + END FUNCTION fomp_in_parallel SUBROUTINE fomp_set_dynamic(dynamic_threads) + USE omp_lib + + IMPLICIT NONE !f2py threadsafe LOGICAL, INTENT(IN) :: dynamic_threads - IF (.FALSE.) PRINT *, dynamic_threads + + + CALL omp_set_dynamic(dynamic_threads) + + + END SUBROUTINE fomp_set_dynamic FUNCTION fomp_get_dynamic() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe LOGICAL :: fomp_get_dynamic - fomp_get_dynamic = .FALSE. + + fomp_get_dynamic = omp_get_dynamic() + + + END FUNCTION fomp_get_dynamic SUBROUTINE fomp_set_nested(nested) + USE omp_lib + + IMPLICIT NONE !f2py threadsafe LOGICAL, INTENT(IN) :: nested - IF (.FALSE.) PRINT *, nested + + + CALL omp_set_nested(nested) + + + END SUBROUTINE fomp_set_nested FUNCTION fomp_get_nested() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe LOGICAL :: fomp_get_nested - fomp_get_nested = .FALSE. + + fomp_get_nested = omp_get_nested() + + + END FUNCTION fomp_get_nested SUBROUTINE fomp_set_schedule(kind, modifier) + USE omp_lib + USE omp_constants, ONLY : fomp_sched_kind IMPLICIT NONE @@ -160,13 +249,20 @@ SUBROUTINE fomp_set_schedule(kind, modifier) INTEGER(KIND=fomp_sched_kind), INTENT(IN) :: kind INTEGER, INTENT(IN) :: modifier - IF (.FALSE.) PRINT *, kind, modifier + + + CALL omp_set_schedule(kind, modifier) + + + END SUBROUTINE fomp_set_schedule SUBROUTINE fomp_get_schedule(kind, modifier) + USE omp_lib + USE omp_constants, ONLY : fomp_sched_kind IMPLICIT NONE @@ -176,121 +272,184 @@ SUBROUTINE fomp_get_schedule(kind, modifier) INTEGER(KIND=fomp_sched_kind), INTENT(OUT) :: kind INTEGER, INTENT(OUT) :: modifier - kind = -1 - modifier = -1 + + CALL omp_get_schedule(kind, modifier) + + + + END SUBROUTINE fomp_get_schedule FUNCTION fomp_get_thread_limit() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER :: fomp_get_thread_limit - fomp_get_thread_limit = -1 + + fomp_get_thread_limit = omp_get_thread_limit() + + + END FUNCTION fomp_get_thread_limit SUBROUTINE fomp_set_max_active_levels(max_levels) + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER, INTENT(IN) :: max_levels - IF (.FALSE.) PRINT *, max_levels + + + CALL omp_set_max_active_levels(max_levels) + + + END SUBROUTINE fomp_set_max_active_levels FUNCTION fomp_get_max_active_levels() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER :: fomp_get_max_active_levels - fomp_get_max_active_levels = -1 + + fomp_get_max_active_levels = omp_get_max_active_levels() + + + END FUNCTION fomp_get_max_active_levels FUNCTION fomp_get_level() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER :: fomp_get_level - fomp_get_level = -1 + + fomp_get_level = omp_get_level() + + + END FUNCTION fomp_get_level FUNCTION fomp_get_ancestor_thread_num(level) + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER, INTENT(IN) :: level INTEGER :: fomp_get_ancestor_thread_num - IF (.FALSE.) PRINT *, level - fomp_get_ancestor_thread_num = -1 + + fomp_get_ancestor_thread_num = omp_get_ancestor_thread_num(level) + + + + END FUNCTION fomp_get_ancestor_thread_num FUNCTION fomp_get_team_size(level) + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER, INTENT(IN) :: level INTEGER :: fomp_get_team_size - IF (.FALSE.) PRINT *, level - fomp_get_team_size = -1 + + fomp_get_team_size = omp_get_team_size(level) + + + + END FUNCTION fomp_get_team_size FUNCTION fomp_get_active_level() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe INTEGER :: fomp_get_active_level - fomp_get_active_level = -1 + + fomp_get_active_level = omp_get_active_level() + + + END FUNCTION fomp_get_active_level FUNCTION fomp_in_final() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe LOGICAL :: fomp_in_final - fomp_in_final = .FALSE. + + fomp_in_final = omp_in_final() + + + END FUNCTION fomp_in_final SUBROUTINE fomp_init_lock(svar) + USE omp_lib + USE omp_constants, ONLY : fomp_lock_kind IMPLICIT NONE @@ -299,13 +458,19 @@ SUBROUTINE fomp_init_lock(svar) INTEGER(KIND=fomp_lock_kind), INTENT(OUT) :: svar - svar = -1 + + CALL omp_init_lock(svar) + + + END SUBROUTINE fomp_init_lock SUBROUTINE fomp_init_nest_lock(nvar) + USE omp_lib + USE omp_constants, ONLY : fomp_nest_lock_kind IMPLICIT NONE @@ -314,13 +479,19 @@ SUBROUTINE fomp_init_nest_lock(nvar) INTEGER(KIND=fomp_nest_lock_kind), INTENT(OUT) :: nvar - nvar = -1 + + CALL omp_init_nest_lock(nvar) + + + END SUBROUTINE fomp_init_nest_lock SUBROUTINE fomp_destroy_lock(svar) + USE omp_lib + USE omp_constants, ONLY : fomp_lock_kind IMPLICIT NONE @@ -328,13 +499,21 @@ SUBROUTINE fomp_destroy_lock(svar) !f2py threadsafe INTEGER(KIND=fomp_lock_kind), INTENT(INOUT) :: svar - IF (.FALSE.) PRINT *, svar + + + CALL omp_destroy_lock(svar) + + + + END SUBROUTINE fomp_destroy_lock SUBROUTINE fomp_destroy_nest_lock(nvar) + USE omp_lib + USE omp_constants, ONLY : fomp_nest_lock_kind IMPLICIT NONE @@ -342,13 +521,20 @@ SUBROUTINE fomp_destroy_nest_lock(nvar) !f2py threadsafe INTEGER(KIND=fomp_nest_lock_kind), INTENT(INOUT) :: nvar - IF (.FALSE.) PRINT *, nvar + + + CALL omp_destroy_nest_lock(nvar) + + + END SUBROUTINE fomp_destroy_nest_lock SUBROUTINE fomp_set_lock(svar) + USE omp_lib + USE omp_constants, ONLY : fomp_lock_kind IMPLICIT NONE @@ -356,13 +542,20 @@ SUBROUTINE fomp_set_lock(svar) !f2py threadsafe INTEGER(KIND=fomp_lock_kind), INTENT(INOUT) :: svar - IF (.FALSE.) PRINT *, svar + + + CALL omp_set_lock(svar) + + + END SUBROUTINE fomp_set_lock SUBROUTINE fomp_set_nest_lock(nvar) + USE omp_lib + USE omp_constants, ONLY : fomp_nest_lock_kind IMPLICIT NONE @@ -370,13 +563,20 @@ SUBROUTINE fomp_set_nest_lock(nvar) !f2py threadsafe INTEGER(KIND=fomp_nest_lock_kind), INTENT(INOUT) :: nvar - IF (.FALSE.) PRINT *, nvar + + + CALL omp_set_nest_lock(nvar) + + + END SUBROUTINE fomp_set_nest_lock SUBROUTINE fomp_unset_lock(svar) + USE omp_lib + USE omp_constants, ONLY : fomp_lock_kind IMPLICIT NONE @@ -384,13 +584,20 @@ SUBROUTINE fomp_unset_lock(svar) !f2py threadsafe INTEGER(KIND=fomp_lock_kind), INTENT(INOUT) :: svar - IF (.FALSE.) PRINT *, svar + + + CALL omp_unset_lock(svar) + + + END SUBROUTINE fomp_unset_lock SUBROUTINE fomp_unset_nest_lock(nvar) + USE omp_lib + USE omp_constants, ONLY : fomp_nest_lock_kind IMPLICIT NONE @@ -398,13 +605,20 @@ SUBROUTINE fomp_unset_nest_lock(nvar) !f2py threadsafe INTEGER(KIND=fomp_nest_lock_kind), INTENT(INOUT) :: nvar - IF (.FALSE.) PRINT *, nvar + + + CALL omp_unset_nest_lock(nvar) + + + END SUBROUTINE fomp_unset_nest_lock FUNCTION fomp_test_lock(svar) + USE omp_lib + USE omp_constants, ONLY : fomp_lock_kind IMPLICIT NONE @@ -413,15 +627,23 @@ FUNCTION fomp_test_lock(svar) INTEGER(KIND=fomp_lock_kind), INTENT(INOUT) :: svar LOGICAL :: fomp_test_lock - IF (.FALSE.) PRINT *, svar - fomp_test_lock = .FALSE. + + fomp_test_lock = omp_test_lock(svar) + + + + + + END FUNCTION fomp_test_lock FUNCTION fomp_test_nest_lock(nvar) + USE omp_lib + USE omp_constants, ONLY : fomp_nest_lock_kind IMPLICIT NONE @@ -430,35 +652,59 @@ FUNCTION fomp_test_nest_lock(nvar) INTEGER(KIND=fomp_nest_lock_kind), INTENT(INOUT) :: nvar INTEGER :: fomp_test_nest_lock - IF (.FALSE.) PRINT *, nvar - fomp_test_nest_lock = -1 + + fomp_test_nest_lock = omp_test_nest_lock(nvar) + + + + + + END FUNCTION fomp_test_nest_lock FUNCTION fomp_get_wtime() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe REAL (KIND=8) :: fomp_get_wtime - fomp_get_wtime = -1 + + fomp_get_wtime = omp_get_wtime() + + + + + END FUNCTION fomp_get_wtime FUNCTION fomp_get_wtick() + USE omp_lib + + IMPLICIT NONE !f2py threadsafe REAL (KIND=8) :: fomp_get_wtick - fomp_get_wtick = -1 + + fomp_get_wtick = omp_get_wtick() + + + + + END FUNCTION fomp_get_wtick diff --git a/fortran/ompgen.F90 b/fortran/ompgen.F90 index 14fa7d4..2ba08b6 100644 --- a/fortran/ompgen.F90 +++ b/fortran/ompgen.F90 @@ -3,7 +3,7 @@ MODULE omp_constants USE omp_lib INTEGER, PARAMETER :: fomp_sched_kind = 4 - INTEGER, PARAMETER :: fomp_lock_kind = 8 + INTEGER, PARAMETER :: fomp_lock_kind = 4 INTEGER, PARAMETER :: fomp_nest_lock_kind = 8 INTEGER(KIND=4), PARAMETER :: fomp_sched_static = 1 INTEGER(KIND=4), PARAMETER :: fomp_sched_dynamic = 2 diff --git a/wrf_ghg_test.ipynb b/wrf_ghg_test.ipynb new file mode 100644 index 0000000..56d098b --- /dev/null +++ b/wrf_ghg_test.ipynb @@ -0,0 +1,674 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + ":241: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility. Expected 16 from C header, got 96 from PyObject\n" + ] + } + ], + "source": [ + "from wrf import getvar\n", + "from netCDF4 import Dataset\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "wrfin = Dataset('./wrfout_d02_2023-07-21_18:00:00')" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "def _extract_ghg(wrffile: Dataset, chem: str):\n", + " if chem == 'xch4':\n", + " _ant = wrffile['CH4_ANT'][0, ...]\n", + " _bck = wrffile['CH4_BCK'][0, ...]\n", + " _bio = wrffile['CH4_BIO'][0, ...]\n", + " elif chem == 'xco2':\n", + " _ant = wrffile['CO2_ANT'][0, ...]\n", + " _bck = wrffile['CO2_BCK'][0, ...]\n", + " _bio = wrffile['CO2_BIO'][0, ...]\n", + " elif chem == 'xco':\n", + " _ant = wrffile['CO_ANT'][0, ...]\n", + " _bck = wrffile['CO_BCK'][0, ...]\n", + " #_bio = wrffile['CO_BIO'][0, ...] ##?\n", + " _bio = np.zeros_like(_bck)\n", + " _ghg = _bio + _ant -_bck\n", + " p = getvar(wrfin, 'P',meta=False)\n", + " sfc_pres = getvar(wrfin, 'PSFC', meta=False)\n", + " if len(_ghg) == len(p):\n", + " pres_bound = np.empty_like(p)\n", + " for i, pres in enumerate(p):\n", + " if i == 0:\n", + " pres_bound[i] = sfc_pres\n", + " pres_bound[i+1] = pres_bound[i] + (2*(pres-pres_bound[i]))\n", + " else:\n", + " try:\n", + " pres_bound[i+1] = pres_bound[i] + (2*(pres-pres_bound[i]))\n", + " except IndexError:\n", + " pass\n", + " p_layer_diff = np.array([pres_bound[i]-pres_bound[i-1] for i in range(1,len(pres_bound))]) #! This may need a value at beginning for xch4[0]\n", + " p_layer_diff = np.insert(p_layer_diff, [0], p[0,...]-sfc_pres, axis=0)\n", + " p_diff = pres_bound[0] - pres_bound[-1]\n", + " return np.sum(_ghg*p_layer_diff)/p_diff" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mainly python: 0.08050990104675293 s\n" + ] + }, + { + "ename": "ValueError", + "evalue": "invalid shape for argument 1 (got (47, 255, 300), expected (255, 300))", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[8], line 7\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mMainly python: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mend\u001b[38;5;241m-\u001b[39mstart\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m s\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 6\u001b[0m start \u001b[38;5;241m=\u001b[39m time\u001b[38;5;241m.\u001b[39mtime()\n\u001b[0;32m----> 7\u001b[0m ghg_getvar \u001b[38;5;241m=\u001b[39m \u001b[43mgetvar\u001b[49m\u001b[43m(\u001b[49m\u001b[43mwrfin\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mxghg\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 8\u001b[0m end \u001b[38;5;241m=\u001b[39m time\u001b[38;5;241m.\u001b[39mtime()\n\u001b[1;32m 9\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgetvar: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mend\u001b[38;5;241m-\u001b[39mstart\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m s\u001b[39m\u001b[38;5;124m'\u001b[39m)\n", + "File \u001b[0;32m~/miniconda3/envs/wrf_python_build/lib/python3.11/site-packages/wrf/routines.py:361\u001b[0m, in \u001b[0;36mgetvar\u001b[0;34m(wrfin, varname, timeidx, method, squeeze, cache, meta, **kwargs)\u001b[0m\n\u001b[1;32m 357\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m is not a valid variable name\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(varname))\n\u001b[1;32m 359\u001b[0m _check_kargs(actual_var, kwargs)\n\u001b[0;32m--> 361\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_FUNC_MAP\u001b[49m\u001b[43m[\u001b[49m\u001b[43mactual_var\u001b[49m\u001b[43m]\u001b[49m\u001b[43m(\u001b[49m\u001b[43mwrfin\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtimeidx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmethod\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msqueeze\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcache\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 362\u001b[0m \u001b[43m \u001b[49m\u001b[43mmeta\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m_key\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/miniconda3/envs/wrf_python_build/lib/python3.11/site-packages/wrf/metadecorators.py:129\u001b[0m, in \u001b[0;36mcopy_and_set_metadata..func_wrapper\u001b[0;34m(wrapped, instance, args, kwargs)\u001b[0m\n\u001b[1;32m 126\u001b[0m new_args, cache_argloc \u001b[38;5;241m=\u001b[39m arg_location(wrapped, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcache\u001b[39m\u001b[38;5;124m\"\u001b[39m, args, kwargs)\n\u001b[1;32m 127\u001b[0m new_args[cache_argloc] \u001b[38;5;241m=\u001b[39m new_cache\n\u001b[0;32m--> 129\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mwrapped\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mnew_args\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 131\u001b[0m outname \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 132\u001b[0m outdimnames \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m()\n", + "File \u001b[0;32m~/miniconda3/envs/wrf_python_build/lib/python3.11/site-packages/wrf/g_ghg.py:100\u001b[0m, in \u001b[0;36mget_xghg\u001b[0;34m(wrfin, timeidx, method, squeeze, cache, meta, _key, ghg)\u001b[0m\n\u001b[1;32m 97\u001b[0m \u001b[38;5;28;01mcase\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01m_\u001b[39;00m:\n\u001b[1;32m 98\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mNotImplementedError\u001b[39;00m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mChemistry variable not supported at this time.\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m--> 100\u001b[0m xghg \u001b[38;5;241m=\u001b[39m \u001b[43m_xghg\u001b[49m\u001b[43m(\u001b[49m\u001b[43msfc_p\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpres\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mant\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbio\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbck\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 102\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m xghg\n", + "File \u001b[0;32m~/miniconda3/envs/wrf_python_build/lib/python3.11/site-packages/wrf/decorators.py:506\u001b[0m, in \u001b[0;36mcheck_args..func_wrapper\u001b[0;34m(wrapped, instance, args, kwargs)\u001b[0m\n\u001b[1;32m 502\u001b[0m \u001b[38;5;66;03m# Check that right dimensions are lined up\u001b[39;00m\n\u001b[1;32m 503\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m (var\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m-\u001b[39mright_var_ndims:] \u001b[38;5;241m!=\u001b[39m\n\u001b[1;32m 504\u001b[0m ref_right_sizes[\u001b[38;5;241m-\u001b[39mright_var_ndims:]):\n\u001b[0;32m--> 506\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minvalid shape for argument \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 507\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m (got \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m, expected \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m)\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(\n\u001b[1;32m 508\u001b[0m i,\n\u001b[1;32m 509\u001b[0m var\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m-\u001b[39mright_var_ndims:],\n\u001b[1;32m 510\u001b[0m ref_right_sizes[\u001b[38;5;241m-\u001b[39mright_var_ndims:]))\n\u001b[1;32m 512\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m wrapped(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n", + "\u001b[0;31mValueError\u001b[0m: invalid shape for argument 1 (got (47, 255, 300), expected (255, 300))" + ] + } + ], + "source": [ + "start = time.time()\n", + "ghg_python = _extract_ghg(wrfin, 'xch4')\n", + "end = time.time()\n", + "print(f'Mainly python: {(end-start):0.4f} s')\n", + "\n", + "start = time.time()\n", + "ghg_getvar = getvar(wrfin, 'xghg')\n", + "end = time.time()\n", + "print(f'getvar: {(end-start):0.4f} s')" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'PSFC' (south_north: 255, west_east: 300)> Size: 306kB\n",
+       "array([[ 99340.37 ,  99345.36 ,  99349.98 , ..., 100734.42 , 100734.35 ,\n",
+       "        100734.67 ],\n",
+       "       [ 99303.21 ,  99311.266,  99315.445, ..., 100732.34 , 100733.12 ,\n",
+       "        100734.19 ],\n",
+       "       [ 99265.07 ,  99274.15 ,  99278.93 , ..., 100732.07 , 100731.44 ,\n",
+       "        100733.836],\n",
+       "       ...,\n",
+       "       [ 93755.04 ,  93721.805,  93690.1  , ...,  99931.79 ,  99947.99 ,\n",
+       "         99968.06 ],\n",
+       "       [ 93765.41 ,  93733.96 ,  93702.164, ...,  99925.1  ,  99941.24 ,\n",
+       "         99960.58 ],\n",
+       "       [ 93777.36 ,  93743.29 ,  93706.53 , ...,  99912.97 ,  99934.305,\n",
+       "         99954.26 ]], dtype=float32)\n",
+       "Coordinates:\n",
+       "    XLONG    (south_north, west_east) float32 306kB -75.46 -75.45 ... -72.07\n",
+       "    XLAT     (south_north, west_east) float32 306kB 40.29 40.29 ... 41.52 41.52\n",
+       "    XTIME    float32 4B 2.52e+03\n",
+       "    Time     datetime64[ns] 8B 2023-07-21T18:00:00\n",
+       "Dimensions without coordinates: south_north, west_east\n",
+       "Attributes:\n",
+       "    FieldType:    104\n",
+       "    MemoryOrder:  XY \n",
+       "    description:  SFC PRESSURE\n",
+       "    units:        Pa\n",
+       "    stagger:      \n",
+       "    coordinates:  XLONG XLAT XTIME\n",
+       "    projection:   LambertConformal(stand_lon=-97.0, moad_cen_lat=40.000007629...
" + ], + "text/plain": [ + " Size: 306kB\n", + "array([[ 99340.37 , 99345.36 , 99349.98 , ..., 100734.42 , 100734.35 ,\n", + " 100734.67 ],\n", + " [ 99303.21 , 99311.266, 99315.445, ..., 100732.34 , 100733.12 ,\n", + " 100734.19 ],\n", + " [ 99265.07 , 99274.15 , 99278.93 , ..., 100732.07 , 100731.44 ,\n", + " 100733.836],\n", + " ...,\n", + " [ 93755.04 , 93721.805, 93690.1 , ..., 99931.79 , 99947.99 ,\n", + " 99968.06 ],\n", + " [ 93765.41 , 93733.96 , 93702.164, ..., 99925.1 , 99941.24 ,\n", + " 99960.58 ],\n", + " [ 93777.36 , 93743.29 , 93706.53 , ..., 99912.97 , 99934.305,\n", + " 99954.26 ]], dtype=float32)\n", + "Coordinates:\n", + " XLONG (south_north, west_east) float32 306kB -75.46 -75.45 ... -72.07\n", + " XLAT (south_north, west_east) float32 306kB 40.29 40.29 ... 41.52 41.52\n", + " XTIME float32 4B 2.52e+03\n", + " Time datetime64[ns] 8B 2023-07-21T18:00:00\n", + "Dimensions without coordinates: south_north, west_east\n", + "Attributes:\n", + " FieldType: 104\n", + " MemoryOrder: XY \n", + " description: SFC PRESSURE\n", + " units: Pa\n", + " stagger: \n", + " coordinates: XLONG XLAT XTIME\n", + " projection: LambertConformal(stand_lon=-97.0, moad_cen_lat=40.000007629..." + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sfc_p = getvar(wrfin, 'PSFC')\n", + "\n", + "sfc_p" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "_wrffortran._wrffortran.dcomputexghg: failed to create array from the 2nd argument `pres` -- 0-th dimension must be fixed to 255 but got 47", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[6], line 13\u001b[0m\n\u001b[1;32m 10\u001b[0m bck \u001b[38;5;241m=\u001b[39m getvar(wrfin, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mCH4_BCK\u001b[39m\u001b[38;5;124m'\u001b[39m, meta\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[1;32m 11\u001b[0m out \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mempty_like(to_np(sfc_p))\n\u001b[0;32m---> 13\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mdcomputexghg\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masfortranarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mto_np\u001b[49m\u001b[43m(\u001b[49m\u001b[43msfc_p\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masfortranarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mp\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masfortranarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mant\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masfortranarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mbio\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masfortranarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mbck\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masfortranarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mout\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 15\u001b[0m \u001b[38;5;28mprint\u001b[39m(result\u001b[38;5;241m.\u001b[39mshape)\n", + "\u001b[0;31mValueError\u001b[0m: _wrffortran._wrffortran.dcomputexghg: failed to create array from the 2nd argument `pres` -- 0-th dimension must be fixed to 255 but got 47" + ] + } + ], + "source": [ + "from wrf._wrffortran import dcomputexghg\n", + "from wrf import to_np, getvar\n", + "import numpy as np\n", + "from netCDF4 import Dataset\n", + "wrfin = Dataset('./wrfout_d02_2023-07-21_18:00:00')\n", + "sfc_p = getvar(wrfin, 'PSFC')\n", + "p = getvar(wrfin, 'P', meta=False)\n", + "ant = getvar(wrfin, 'CH4_ANT', meta=False)\n", + "bio = getvar(wrfin, 'CH4_BIO', meta=False)\n", + "bck = getvar(wrfin, 'CH4_BCK', meta=False)\n", + "out = np.empty_like(to_np(sfc_p))\n", + "\n", + "result = dcomputexghg(np.asfortranarray(to_np(sfc_p)), np.asfortranarray(p), np.asfortranarray(ant), np.asfortranarray(bio), np.asfortranarray(bck), np.asfortranarray(out))\n", + "\n", + "print(result.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(47, 255, 300)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "p.shape" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "wrf_python_build", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.10" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}