diff --git a/.github/workflows/integration-tests-ci.yml b/.github/workflows/integration-tests-ci.yml new file mode 100644 index 0000000000..c841352b0d --- /dev/null +++ b/.github/workflows/integration-tests-ci.yml @@ -0,0 +1,60 @@ +name: Integration Tests + +# The CloudSC corpus tests parse a multi-thousand-block SDFG and compile it, which is far too +# slow for the per-push ``General Tests`` matrix -- they carry the ``long`` marker that matrix +# excludes, so without this workflow nothing runs them. Kept to a single Python version: the +# job is bound by the parse and the C++ compile, not by anything version-specific. + +on: + push: + branches: [ main, ci-fix ] + pull_request: + branches: [ main, ci-fix ] + merge_group: + branches: [ main, ci-fix ] + workflow_dispatch: + +concurrency: + group: ${{github.workflow}}-${{github.ref}} + cancel-in-progress: true + +jobs: + cloudsc-e2e: + if: "!contains(github.event.pull_request.labels.*.name, 'no-ci')" + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ['3.14'] + + steps: + - uses: actions/checkout@v7 + with: + submodules: 'recursive' + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v6 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + # Make dependency setup faster + echo 'set man-db/auto-update false' | sudo debconf-communicate >/dev/null + sudo dpkg-reconfigure man-db + # Install dependencies + sudo apt-get update + sudo apt-get install -y libyaml-dev cmake libblas-dev libopenblas-dev liblapacke-dev + pip install -e ".[testing]" + + - name: CloudSC corpus tests + run: | + export NOSTATUSBAR=1 + export DACE_cache=unique + # The tests build with ``simplify=False`` and apply simplification as an explicit step, + # so the automatic heuristic must not have run first. + export DACE_optimizer_automatic_simplification=0 + # Ask the installed package where it lives instead of assuming the job's directory. + ROOT=$(python -c 'import pathlib, dace; print(pathlib.Path(dace.__file__).resolve().parents[1])') + # No xdist: the SDFG is parsed once per process and memoized, so a worker per test would + # pay the multi-minute parse again instead of sharing one. + pytest --tb=short --timeout_method thread --timeout=3600 \ + "$ROOT/tests/corpus/cloudsc_regression_test.py" \ + "$ROOT/tests/passes/constant_propagation_on_cloudsc_test.py" diff --git a/tests/corpus/__init__.py b/tests/corpus/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/corpus/cloudsc/README.md b/tests/corpus/cloudsc/README.md new file mode 100644 index 0000000000..ab72854728 --- /dev/null +++ b/tests/corpus/cloudsc/README.md @@ -0,0 +1,24 @@ +# CloudSC test corpus + +The ECMWF `dwarf-p-cloudsc` cloud microphysics kernel, inlined into a single `dace.program` +(`cloudsc.py`). It is callback-free, so `cloudsc_py.to_sdfg()` builds standalone. The result is a +large, wide, deeply nested SDFG (thousands of blocks, nested loop regions many levels deep), which +makes it useful as a scaling test for whole-SDFG analyses and passes. + +`generate_data_for_cloudsc.py` provides: + +- `build_cloudsc_sdfg(simplify=False)` — the parsed SDFG. The parse takes minutes, so it is memoized + per process in `PARSED_CLOUDSC` and every caller gets a deepcopy: an SDFG is mutable and every + consumer transforms it, so handing out the memoized object would leak one test's edits into the + next. Nothing is written to disk, and under pytest-xdist each worker pays the parse once. +- `generate_cloudsc_inputs(sdfg, seed)` — a physically realistic input set. The physical constants + and the per-array `[min, max]` ranges are the values from the dwarf's `config-files/input.h5`, + mirrored here so nothing external is needed. Random inputs would sit on every threshold in the + kernel, where a harmless floating-point reassociation flips a branch and looks like a bug. +- `run_and_compare(reference, candidate)` — drives both SDFGs on the same inputs and compares every + output array. With the default IEEE build (`-O0`, no fast-math, no FP contraction) and sequential + schedules, a value-preserving transformation reproduces the reference bit-for-bit, hence the + `1e-15` default tolerance. + +The grid is small (`klev = klon = 32`) so a compiled run is quick; the input ranges are bounds, not +vertical profiles, so they stay valid at any size. diff --git a/tests/corpus/cloudsc/__init__.py b/tests/corpus/cloudsc/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/corpus/cloudsc/cloudsc.py b/tests/corpus/cloudsc/cloudsc.py new file mode 100644 index 0000000000..e30408e14d --- /dev/null +++ b/tests/corpus/cloudsc/cloudsc.py @@ -0,0 +1,1363 @@ +# Copyright 2019-2026 ETH Zurich and the DaCe authors. All rights reserved. +"""Inlined CloudSC (ECMWF dwarf-p-cloudsc) cloud microphysics kernel as a single ``dace.program``. + +Callback-free and self-contained, so ``cloudsc_py.to_sdfg()`` yields a large, wide, deeply nested +SDFG that compiles and runs standalone. Input data generation lives in +``generate_data_for_cloudsc.py``. +""" +import numpy as np +import dace + +klon = dace.symbol('klon', dtype=dace.int32) +klev = dace.symbol('klev', dtype=dace.int32) +nclv = dace.symbol('nclv', dtype=dace.int32) +ncldql = dace.symbol('ncldql', dtype=dace.int32) +ncldqi = dace.symbol('ncldqi', dtype=dace.int32) +ncldqr = dace.symbol('ncldqr', dtype=dace.int32) +ncldqs = dace.symbol('ncldqs', dtype=dace.int32) +ncldqv = dace.symbol('ncldqv', dtype=dace.int32) + + +@dace.program +def cloudsc_py( + kidia: dace.int32, + kfdia: dace.int32, + ptsphy: dace.float64, + pt: dace.float64[klev, klon], + pq: dace.float64[klev, klon], + tendency_tmp_t: dace.float64[klev, klon], + tendency_tmp_q: dace.float64[klev, klon], + tendency_tmp_a: dace.float64[klev, klon], + tendency_tmp_cld: dace.float64[nclv, klev, klon], + tendency_loc_t: dace.float64[klev, klon], + tendency_loc_q: dace.float64[klev, klon], + tendency_loc_a: dace.float64[klev, klon], + tendency_loc_cld: dace.float64[nclv, klev, klon], + pvfa: dace.float64[klev, klon], + pvfl: dace.float64[klev, klon], + pvfi: dace.float64[klev, klon], + pdyna: dace.float64[klev, klon], + pdynl: dace.float64[klev, klon], + pdyni: dace.float64[klev, klon], + phrsw: dace.float64[klev, klon], + phrlw: dace.float64[klev, klon], + pvervel: dace.float64[klev, klon], + pap: dace.float64[klev, klon], + paph: dace.float64[klev + 1, klon], + plsm: dace.float64[klon], + ldcum: dace.int32[klon], + ktype: dace.int32[klon], + plu: dace.float64[klev, klon], + plude: dace.float64[klev, klon], + psnde: dace.float64[klev, klon], + pmfu: dace.float64[klev, klon], + pmfd: dace.float64[klev, klon], + pa: dace.float64[klev, klon], + pclv: dace.float64[nclv, klev, klon], + psupsat: dace.float64[klev, klon], + plcrit_aer: dace.float64[klev, klon], + picrit_aer: dace.float64[klev, klon], + pre_ice: dace.float64[klev, klon], + pccn: dace.float64[klev, klon], + pnice: dace.float64[klev, klon], + pcovptot: dace.float64[klev, klon], + prainfrac_toprfz: dace.float64[klon], + pfsqlf: dace.float64[klev + 1, klon], + pfsqif: dace.float64[klev + 1, klon], + pfcqnng: dace.float64[klev + 1, klon], + pfcqlng: dace.float64[klev + 1, klon], + pfsqrf: dace.float64[klev + 1, klon], + pfsqsf: dace.float64[klev + 1, klon], + pfcqrng: dace.float64[klev + 1, klon], + pfcqsng: dace.float64[klev + 1, klon], + pfsqltur: dace.float64[klev + 1, klon], + pfsqitur: dace.float64[klev + 1, klon], + pfplsl: dace.float64[klev + 1, klon], + pfplsn: dace.float64[klev + 1, klon], + pfhpsl: dace.float64[klev + 1, klon], + pfhpsn: dace.float64[klev + 1, klon], + # --- YDCST (flattened) --- + ydcst_rg: dace.float64, + ydcst_rd: dace.float64, + ydcst_rcpd: dace.float64, + ydcst_retv: dace.float64, + ydcst_rlvtt: dace.float64, + ydcst_rlstt: dace.float64, + ydcst_rlmlt: dace.float64, + ydcst_rtt: dace.float64, + ydcst_rv: dace.float64, + # --- YDTHF (flattened) --- + ydthf_r2es: dace.float64, + ydthf_r3les: dace.float64, + ydthf_r3ies: dace.float64, + ydthf_r4les: dace.float64, + ydthf_r4ies: dace.float64, + ydthf_r5les: dace.float64, + ydthf_r5ies: dace.float64, + ydthf_r5alvcp: dace.float64, + ydthf_r5alscp: dace.float64, + ydthf_ralvdcp: dace.float64, + ydthf_ralsdcp: dace.float64, + ydthf_ralfdcp: dace.float64, + ydthf_rtwat: dace.float64, + ydthf_rtice: dace.float64, + ydthf_rticecu: dace.float64, + ydthf_rtwat_rtice_r: dace.float64, + ydthf_rtwat_rticecu_r: dace.float64, + ydthf_rkoop1: dace.float64, + ydthf_rkoop2: dace.float64, + # --- YRECLDP (flattened) --- + yrecldp_ramid: dace.float64, + yrecldp_rcldiff: dace.float64, + yrecldp_rcldiff_convi: dace.float64, + yrecldp_ramin: dace.float64, + yrecldp_rlmin: dace.float64, + yrecldp_rdensref: dace.float64, + yrecldp_rtaumel: dace.float64, + yrecldp_rvice: dace.float64, + yrecldp_rvrain: dace.float64, + yrecldp_rvsnow: dace.float64, + yrecldp_rthomo: dace.float64, + yrecldp_rcovpmin: dace.float64, + yrecldp_rkooptau: dace.float64, + yrecldp_rcldtopcf: dace.float64, + yrecldp_rkconv: dace.float64, + yrecldp_rclcrit_land: dace.float64, + yrecldp_rclcrit_sea: dace.float64, + yrecldp_rlcritsnow: dace.float64, + yrecldp_rprecrhmax: dace.float64, + yrecldp_rprc1: dace.float64, + yrecldp_rvrfactor: dace.float64, + yrecldp_rpecons: dace.float64, + yrecldp_rnice: dace.float64, + yrecldp_riceinit: dace.float64, + yrecldp_rdepliqrefrate: dace.float64, + yrecldp_rdepliqrefdepth: dace.float64, + yrecldp_rsnowlin1: dace.float64, + yrecldp_rsnowlin2: dace.float64, + yrecldp_rccn: dace.float64, + yrecldp_nssopt: dace.int32, + yrecldp_ncldtop: dace.int32, + yrecldp_laericesed: dace.int32, + yrecldp_laerliqautolsp: dace.int32, + yrecldp_laerliqcoll: dace.int32, + yrecldp_laericeauto: dace.int32, + # --- YRECLDP RCL_* microphysics constants --- + yrecldp_rcl_kkaau: dace.float64, + yrecldp_rcl_kkbauq: dace.float64, + yrecldp_rcl_kkbaun: dace.float64, + yrecldp_rcl_kkaac: dace.float64, + yrecldp_rcl_kkbac: dace.float64, + yrecldp_rcl_kk_cloud_num_land: dace.float64, + yrecldp_rcl_kk_cloud_num_sea: dace.float64, + yrecldp_rcl_fac1: dace.float64, + yrecldp_rcl_fac2: dace.float64, + yrecldp_rcl_fzrab: dace.float64, + yrecldp_rcl_apb1: dace.float64, + yrecldp_rcl_apb2: dace.float64, + yrecldp_rcl_apb3: dace.float64, + yrecldp_rcl_const1i: dace.float64, + yrecldp_rcl_const2i: dace.float64, + yrecldp_rcl_const3i: dace.float64, + yrecldp_rcl_const4i: dace.float64, + yrecldp_rcl_const5i: dace.float64, + yrecldp_rcl_const6i: dace.float64, + yrecldp_rcl_const1s: dace.float64, + yrecldp_rcl_const2s: dace.float64, + yrecldp_rcl_const3s: dace.float64, + yrecldp_rcl_const4s: dace.float64, + yrecldp_rcl_const5s: dace.float64, + yrecldp_rcl_const6s: dace.float64, + yrecldp_rcl_const7s: dace.float64, + yrecldp_rcl_const8s: dace.float64, + yrecldp_rcl_const1r: dace.float64, + yrecldp_rcl_const2r: dace.float64, + yrecldp_rcl_const3r: dace.float64, + yrecldp_rcl_const4r: dace.float64, + yrecldp_rcl_const5r: dace.float64, + yrecldp_rcl_const6r: dace.float64, + yrecldp_rcl_ka273: dace.float64, + yrecldp_rcl_cdenom1: dace.float64, + yrecldp_rcl_cdenom2: dace.float64, + yrecldp_rcl_cdenom3: dace.float64, +): + zlcond1 = np.ndarray(shape=(klon, ), dtype=np.float64) + zlcond2 = np.ndarray(shape=(klon, ), dtype=np.float64) + zlevapl = np.ndarray(shape=(klon, ), dtype=np.float64) + zlevapi = np.ndarray(shape=(klon, ), dtype=np.float64) + zrainaut = np.ndarray(shape=(klon, ), dtype=np.float64) + zsnowaut = np.ndarray(shape=(klon, ), dtype=np.float64) + zliqcld = np.ndarray(shape=(klon, ), dtype=np.float64) + zicecld = np.ndarray(shape=(klon, ), dtype=np.float64) + zfokoop = np.ndarray(shape=(klon, ), dtype=np.float64) + zicenuclei = np.ndarray(shape=(klon, ), dtype=np.float64) + zlicld = np.ndarray(shape=(klon, ), dtype=np.float64) + zlfinalsum = np.ndarray(shape=(klon, ), dtype=np.float64) + zdqs = np.ndarray(shape=(klon, ), dtype=np.float64) + ztold = np.ndarray(shape=(klon, ), dtype=np.float64) + zqold = np.ndarray(shape=(klon, ), dtype=np.float64) + zdtgdp = np.ndarray(shape=(klon, ), dtype=np.float64) + zrdtgdp = np.ndarray(shape=(klon, ), dtype=np.float64) + ztrpaus = np.ndarray(shape=(klon, ), dtype=np.float64) + zcovpclr = np.ndarray(shape=(klon, ), dtype=np.float64) + zcovptot = np.ndarray(shape=(klon, ), dtype=np.float64) + zcovpmax = np.ndarray(shape=(klon, ), dtype=np.float64) + zqpretot = np.ndarray(shape=(klon, ), dtype=np.float64) + zldefr = np.ndarray(shape=(klon, ), dtype=np.float64) + zldifdt = np.ndarray(shape=(klon, ), dtype=np.float64) + zdtgdpf = np.ndarray(shape=(klon, ), dtype=np.float64) + zacust = np.ndarray(shape=(klon, ), dtype=np.float64) + zmf = np.ndarray(shape=(klon, ), dtype=np.float64) + zrho = np.ndarray(shape=(klon, ), dtype=np.float64) + ztmp1 = np.ndarray(shape=(klon, ), dtype=np.float64) + ztmp2 = np.ndarray(shape=(klon, ), dtype=np.float64) + ztmp3 = np.ndarray(shape=(klon, ), dtype=np.float64) + ztmp4 = np.ndarray(shape=(klon, ), dtype=np.float64) + ztmp5 = np.ndarray(shape=(klon, ), dtype=np.float64) + ztmp6 = np.ndarray(shape=(klon, ), dtype=np.float64) + ztmp7 = np.ndarray(shape=(klon, ), dtype=np.float64) + zalfawm = np.ndarray(shape=(klon, ), dtype=np.float64) + zsolab = np.ndarray(shape=(klon, ), dtype=np.float64) + zsolac = np.ndarray(shape=(klon, ), dtype=np.float64) + zanewm1 = np.ndarray(shape=(klon, ), dtype=np.float64) + zgdp = np.ndarray(shape=(klon, ), dtype=np.float64) + zda = np.ndarray(shape=(klon, ), dtype=np.float64) + zdp = np.ndarray(shape=(klon, ), dtype=np.float64) + zpaphd = np.ndarray(shape=(klon, ), dtype=np.float64) + zmin = np.ndarray(shape=(klon, ), dtype=np.float64) + zsupsat = np.ndarray(shape=(klon, ), dtype=np.float64) + zmeltmax = np.ndarray(shape=(klon, ), dtype=np.float64) + zfrzmax = np.ndarray(shape=(klon, ), dtype=np.float64) + zicetot = np.ndarray(shape=(klon, ), dtype=np.float64) + zdqsliqdt = np.ndarray(shape=(klon, ), dtype=np.float64) + zdqsicedt = np.ndarray(shape=(klon, ), dtype=np.float64) + zdqsmixdt = np.ndarray(shape=(klon, ), dtype=np.float64) + zcorqsliq = np.ndarray(shape=(klon, ), dtype=np.float64) + zcorqsice = np.ndarray(shape=(klon, ), dtype=np.float64) + zcorqsmix = np.ndarray(shape=(klon, ), dtype=np.float64) + zevaplimliq = np.ndarray(shape=(klon, ), dtype=np.float64) + zevaplimice = np.ndarray(shape=(klon, ), dtype=np.float64) + zevaplimmix = np.ndarray(shape=(klon, ), dtype=np.float64) + zcldtopdist = np.ndarray(shape=(klon, ), dtype=np.float64) + zrainacc = np.ndarray(shape=(klon, ), dtype=np.float64) + zraincld = np.ndarray(shape=(klon, ), dtype=np.float64) + zsnowrime = np.ndarray(shape=(klon, ), dtype=np.float64) + zsnowcld = np.ndarray(shape=(klon, ), dtype=np.float64) + zrg = np.ndarray(shape=(klon, ), dtype=np.float64) + psum_solqa = np.ndarray(shape=(klon, ), dtype=np.float64) + llflag = np.ndarray(shape=(klon, ), dtype=np.float64) + llrainliq = np.ndarray(shape=(klon, ), dtype=np.int32) + iphase = np.ndarray(shape=(nclv, ), dtype=np.int32) + imelt = np.ndarray(shape=(nclv, ), dtype=np.int32) + llfall = np.ndarray(shape=(nclv, ), dtype=np.int32) + zvqx = np.ndarray(shape=(nclv, ), dtype=np.float64) + zfoealfa = np.ndarray(shape=(klev + 1, klon), dtype=np.float64) + ztp1 = np.ndarray(shape=(klev, klon), dtype=np.float64) + zlcust = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zli = np.ndarray(shape=(klev, klon), dtype=np.float64) + za = np.ndarray(shape=(klev, klon), dtype=np.float64) + zaorig = np.ndarray(shape=(klev, klon), dtype=np.float64) + llindex1 = np.ndarray(shape=(nclv, klon), dtype=np.int32) + llindex3 = np.ndarray(shape=(nclv, nclv, klon), dtype=np.int32) + iorder = np.ndarray(shape=(nclv, klon), dtype=np.int32) + zliqfrac = np.ndarray(shape=(klev, klon), dtype=np.float64) + zicefrac = np.ndarray(shape=(klev, klon), dtype=np.float64) + zqx = np.ndarray(shape=(nclv, klev, klon), dtype=np.float64) + zqx0 = np.ndarray(shape=(nclv, klev, klon), dtype=np.float64) + zqxn = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zqxfg = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zqxnm1 = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zfluxq = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zpfplsx = np.ndarray(shape=(nclv, klev + 1, klon), dtype=np.float64) + zlneg = np.ndarray(shape=(nclv, klev, klon), dtype=np.float64) + zqxn2d = np.ndarray(shape=(nclv, klev, klon), dtype=np.float64) + zqsmix = np.ndarray(shape=(klev, klon), dtype=np.float64) + zqsliq = np.ndarray(shape=(klev, klon), dtype=np.float64) + zqsice = np.ndarray(shape=(klev, klon), dtype=np.float64) + zfoeewmt = np.ndarray(shape=(klev, klon), dtype=np.float64) + zfoeew = np.ndarray(shape=(klev, klon), dtype=np.float64) + zfoeeliqt = np.ndarray(shape=(klev, klon), dtype=np.float64) + zsolqa = np.ndarray(shape=(nclv, nclv, klon), dtype=np.float64) + zsolqb = np.ndarray(shape=(nclv, nclv, klon), dtype=np.float64) + zqlhs = np.ndarray(shape=(nclv, nclv, klon), dtype=np.float64) + zratio = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zsinksum = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zfallsink = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zfallsrce = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zconvsrce = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zconvsink = np.ndarray(shape=(nclv, klon), dtype=np.float64) + zpsupsatsrce = np.ndarray(shape=(nclv, klon), dtype=np.float64) + ztw1 = 1329.31 + ztw2 = 0.0074615 + ztw3 = 85000.0 + ztw4 = 40.637 + ztw5 = 275.0 + zepsilon = 1e-14 + iwarmrain = 2 + ievaprain = 2 + ievapsnow = 1 + idepice = 1 + zqtmst = 1.0 / ptsphy + zgdcp = ydcst_rg / ydcst_rcpd + zrdcp = ydcst_rd / ydcst_rcpd + zcons1a = ydcst_rcpd / (ydcst_rlmlt * ydcst_rg * yrecldp_rtaumel) + zepsec = 1e-14 + zrg_r = 1.0 / ydcst_rg + zrldcp = 1.0 / (ydthf_ralsdcp - ydthf_ralvdcp) + iphase[ncldqv - 1] = 0 + iphase[ncldql - 1] = 1 + iphase[ncldqr - 1] = 1 + iphase[ncldqi - 1] = 2 + iphase[ncldqs - 1] = 2 + imelt[ncldqv - 1] = -99 + imelt[ncldql - 1] = ncldqi + imelt[ncldqr - 1] = ncldqs + imelt[ncldqi - 1] = ncldqr + imelt[ncldqs - 1] = ncldqr + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + tendency_loc_t[jk - 1, jl - 1] = 0.0 + tendency_loc_q[jk - 1, jl - 1] = 0.0 + tendency_loc_a[jk - 1, jl - 1] = 0.0 + for jm in range(1, nclv - 1 + 1): + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + tendency_loc_cld[jm - 1, jk - 1, jl - 1] = 0.0 + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + pcovptot[jk - 1, jl - 1] = 0.0 + tendency_loc_cld[nclv - 1, jk - 1, jl - 1] = 0.0 + zvqx[ncldqv - 1] = 0.0 + zvqx[ncldql - 1] = 0.0 + zvqx[ncldqi - 1] = yrecldp_rvice + zvqx[ncldqr - 1] = yrecldp_rvrain + zvqx[ncldqs - 1] = yrecldp_rvsnow + for jm in range(1, nclv + 1): + llfall[jm - 1] = False + for jm in range(1, nclv + 1): + if zvqx[jm - 1] > 0.0: + llfall[jm - 1] = True + llfall[ncldqi - 1] = False + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + ztp1[jk - 1, jl - 1] = pt[jk - 1, jl - 1] + ptsphy * tendency_tmp_t[jk - 1, jl - 1] + zqx[ncldqv - 1, jk - 1, jl - 1] = pq[jk - 1, jl - 1] + ptsphy * tendency_tmp_q[jk - 1, jl - 1] + zqx0[ncldqv - 1, jk - 1, jl - 1] = pq[jk - 1, jl - 1] + ptsphy * tendency_tmp_q[jk - 1, jl - 1] + za[jk - 1, jl - 1] = pa[jk - 1, jl - 1] + ptsphy * tendency_tmp_a[jk - 1, jl - 1] + zaorig[jk - 1, jl - 1] = pa[jk - 1, jl - 1] + ptsphy * tendency_tmp_a[jk - 1, jl - 1] + for jm in range(1, nclv - 1 + 1): + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + zqx[jm - 1, jk - 1, + jl - 1] = pclv[jm - 1, jk - 1, jl - 1] + ptsphy * tendency_tmp_cld[jm - 1, jk - 1, jl - 1] + zqx0[jm - 1, jk - 1, + jl - 1] = pclv[jm - 1, jk - 1, jl - 1] + ptsphy * tendency_tmp_cld[jm - 1, jk - 1, jl - 1] + for jm in range(1, nclv + 1): + for jk in range(1, klev + 1 + 1): + for jl in range(kidia, kfdia + 1): + zpfplsx[jm - 1, jk - 1, jl - 1] = 0.0 + for jm in range(1, nclv + 1): + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + zqxn2d[jm - 1, jk - 1, jl - 1] = 0.0 + zlneg[jm - 1, jk - 1, jl - 1] = 0.0 + for jl in range(kidia, kfdia + 1): + prainfrac_toprfz[jl - 1] = 0.0 + for jl in range(1, klon + 1): + llrainliq[jl - 1] = True + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + if zqx[ncldql - 1, jk - 1, jl - 1] + zqx[ncldqi - 1, jk - 1, + jl - 1] < yrecldp_rlmin or za[jk - 1, jl - 1] < yrecldp_ramin: + zlneg[ncldql - 1, jk - 1, jl - 1] = zlneg[ncldql - 1, jk - 1, jl - 1] + zqx[ncldql - 1, jk - 1, jl - 1] + zqadj = zqx[ncldql - 1, jk - 1, jl - 1] * zqtmst + tendency_loc_q[jk - 1, jl - 1] = tendency_loc_q[jk - 1, jl - 1] + zqadj + tendency_loc_t[jk - 1, jl - 1] = tendency_loc_t[jk - 1, jl - 1] - ydthf_ralvdcp * zqadj + zqx[ncldqv - 1, jk - 1, jl - 1] = zqx[ncldqv - 1, jk - 1, jl - 1] + zqx[ncldql - 1, jk - 1, jl - 1] + zqx[ncldql - 1, jk - 1, jl - 1] = 0.0 + zlneg[ncldqi - 1, jk - 1, jl - 1] = zlneg[ncldqi - 1, jk - 1, jl - 1] + zqx[ncldqi - 1, jk - 1, jl - 1] + zqadj = zqx[ncldqi - 1, jk - 1, jl - 1] * zqtmst + tendency_loc_q[jk - 1, jl - 1] = tendency_loc_q[jk - 1, jl - 1] + zqadj + tendency_loc_t[jk - 1, jl - 1] = tendency_loc_t[jk - 1, jl - 1] - ydthf_ralsdcp * zqadj + zqx[ncldqv - 1, jk - 1, jl - 1] = zqx[ncldqv - 1, jk - 1, jl - 1] + zqx[ncldqi - 1, jk - 1, jl - 1] + zqx[ncldqi - 1, jk - 1, jl - 1] = 0.0 + za[jk - 1, jl - 1] = 0.0 + for jm in range(1, nclv - 1 + 1): + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + if zqx[jm - 1, jk - 1, jl - 1] < yrecldp_rlmin: + zlneg[jm - 1, jk - 1, jl - 1] = zlneg[jm - 1, jk - 1, jl - 1] + zqx[jm - 1, jk - 1, jl - 1] + zqadj = zqx[jm - 1, jk - 1, jl - 1] * zqtmst + tendency_loc_q[jk - 1, jl - 1] = tendency_loc_q[jk - 1, jl - 1] + zqadj + if iphase[jm - 1] == 1: + tendency_loc_t[jk - 1, jl - 1] = tendency_loc_t[jk - 1, jl - 1] - ydthf_ralvdcp * zqadj + if iphase[jm - 1] == 2: + tendency_loc_t[jk - 1, jl - 1] = tendency_loc_t[jk - 1, jl - 1] - ydthf_ralsdcp * zqadj + zqx[ncldqv - 1, jk - 1, jl - 1] = zqx[ncldqv - 1, jk - 1, jl - 1] + zqx[jm - 1, jk - 1, jl - 1] + zqx[jm - 1, jk - 1, jl - 1] = 0.0 + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + zfoealfa[jk - 1, jl - 1] = min(1.0, + ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2) + zfoeewmt[jk - 1, jl - 1] = min( + ydthf_r2es * + (min(1.0, + ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * ydthf_rtwat_rtice_r)** + 2) * np.exp(ydthf_r3les * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4les)) + + (1.0 - min(1.0, ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2)) * np.exp(ydthf_r3ies * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4ies))) / + pap[jk - 1, jl - 1], 0.5) + zqsmix[jk - 1, jl - 1] = zfoeewmt[jk - 1, jl - 1] + zqsmix[jk - 1, jl - 1] = zqsmix[jk - 1, jl - 1] / (1.0 - ydcst_retv * zqsmix[jk - 1, jl - 1]) + zalfa = max(0.0, 1.0 * np.sign(ztp1[jk - 1, jl - 1] - ydcst_rtt)) + zfoeew[jk - 1, jl - 1] = min( + (zalfa * (ydthf_r2es * np.exp(ydthf_r3les * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4les))) + (1.0 - zalfa) * + (ydthf_r2es * np.exp(ydthf_r3ies * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4ies)))) / pap[jk - 1, jl - 1], 0.5) + zfoeew[jk - 1, jl - 1] = min(0.5, zfoeew[jk - 1, jl - 1]) + zqsice[jk - 1, jl - 1] = zfoeew[jk - 1, jl - 1] / (1.0 - ydcst_retv * zfoeew[jk - 1, jl - 1]) + zfoeeliqt[jk - 1, jl - 1] = min( + ydthf_r2es * np.exp(ydthf_r3les * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4les)) / pap[jk - 1, jl - 1], 0.5) + zqsliq[jk - 1, jl - 1] = zfoeeliqt[jk - 1, jl - 1] + zqsliq[jk - 1, jl - 1] = zqsliq[jk - 1, jl - 1] / (1.0 - ydcst_retv * zqsliq[jk - 1, jl - 1]) + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + za[jk - 1, jl - 1] = max(0.0, min(1.0, za[jk - 1, jl - 1])) + zli[jk - 1, jl - 1] = zqx[ncldql - 1, jk - 1, jl - 1] + zqx[ncldqi - 1, jk - 1, jl - 1] + if zli[jk - 1, jl - 1] > yrecldp_rlmin: + zliqfrac[jk - 1, jl - 1] = zqx[ncldql - 1, jk - 1, jl - 1] / zli[jk - 1, jl - 1] + zicefrac[jk - 1, jl - 1] = 1.0 - zliqfrac[jk - 1, jl - 1] + else: + zliqfrac[jk - 1, jl - 1] = 0.0 + zicefrac[jk - 1, jl - 1] = 0.0 + for jl in range(kidia, kfdia + 1): + ztrpaus[jl - 1] = 0.1 + zpaphd[jl - 1] = 1.0 / paph[klev + 1 - 1, jl - 1] + for jk in range(1, klev - 1 + 1): + for jl in range(kidia, kfdia + 1): + zsig = pap[jk - 1, jl - 1] * zpaphd[jl - 1] + if zsig > 0.1 and zsig < 0.4 and (ztp1[jk - 1, jl - 1] > ztp1[jk + 1 - 1, jl - 1]): + ztrpaus[jl - 1] = zsig + for jl in range(kidia, kfdia + 1): + zanewm1[jl - 1] = 0.0 + zda[jl - 1] = 0.0 + zcovpclr[jl - 1] = 0.0 + zcovpmax[jl - 1] = 0.0 + zcovptot[jl - 1] = 0.0 + zcldtopdist[jl - 1] = 0.0 + for jk in range(yrecldp_ncldtop, klev + 1): + for jm in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zqxfg[jm - 1, jl - 1] = zqx[jm - 1, jk - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + zlicld[jl - 1] = 0.0 + zrainaut[jl - 1] = 0.0 + zrainacc[jl - 1] = 0.0 + zsnowaut[jl - 1] = 0.0 + zldefr[jl - 1] = 0.0 + zacust[jl - 1] = 0.0 + zqpretot[jl - 1] = 0.0 + zlfinalsum[jl - 1] = 0.0 + zlcond1[jl - 1] = 0.0 + zlcond2[jl - 1] = 0.0 + zsupsat[jl - 1] = 0.0 + zlevapl[jl - 1] = 0.0 + zlevapi[jl - 1] = 0.0 + zsolab[jl - 1] = 0.0 + zsolac[jl - 1] = 0.0 + zicetot[jl - 1] = 0.0 + for jm in range(1, nclv + 1): + for jn in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zsolqb[jm - 1, jn - 1, jl - 1] = 0.0 + zsolqa[jm - 1, jn - 1, jl - 1] = 0.0 + for jm in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zfallsrce[jm - 1, jl - 1] = 0.0 + zfallsink[jm - 1, jl - 1] = 0.0 + zconvsrce[jm - 1, jl - 1] = 0.0 + zconvsink[jm - 1, jl - 1] = 0.0 + zpsupsatsrce[jm - 1, jl - 1] = 0.0 + zratio[jm - 1, jl - 1] = 0.0 + for jl in range(kidia, kfdia + 1): + zdp[jl - 1] = paph[jk + 1 - 1, jl - 1] - paph[jk - 1, jl - 1] + zgdp[jl - 1] = ydcst_rg / zdp[jl - 1] + zrho[jl - 1] = pap[jk - 1, jl - 1] / (ydcst_rd * ztp1[jk - 1, jl - 1]) + zdtgdp[jl - 1] = ptsphy * zgdp[jl - 1] + zrdtgdp[jl - 1] = zdp[jl - 1] * (1.0 / (ptsphy * ydcst_rg)) + if jk > 1: + zdtgdpf[jl - 1] = ptsphy * ydcst_rg / (pap[jk - 1, jl - 1] - pap[jk - 1 - 1, jl - 1]) + zfacw = ydthf_r5les / (ztp1[jk - 1, jl - 1] - ydthf_r4les)**2 + zcor = 1.0 / (1.0 - ydcst_retv * zfoeeliqt[jk - 1, jl - 1]) + zdqsliqdt[jl - 1] = zfacw * zcor * zqsliq[jk - 1, jl - 1] + zcorqsliq[jl - 1] = 1.0 + ydthf_ralvdcp * zdqsliqdt[jl - 1] + zfaci = ydthf_r5ies / (ztp1[jk - 1, jl - 1] - ydthf_r4ies)**2 + zcor = 1.0 / (1.0 - ydcst_retv * zfoeew[jk - 1, jl - 1]) + zdqsicedt[jl - 1] = zfaci * zcor * zqsice[jk - 1, jl - 1] + zcorqsice[jl - 1] = 1.0 + ydthf_ralsdcp * zdqsicedt[jl - 1] + zalfaw = zfoealfa[jk - 1, jl - 1] + zalfawm[jl - 1] = zalfaw + zfac = zalfaw * zfacw + (1.0 - zalfaw) * zfaci + zcor = 1.0 / (1.0 - ydcst_retv * zfoeewmt[jk - 1, jl - 1]) + zdqsmixdt[jl - 1] = zfac * zcor * zqsmix[jk - 1, jl - 1] + zcorqsmix[jl - + 1] = 1.0 + (min(1.0, ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2) * ydthf_ralvdcp + + (1.0 - min(1.0, + ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2)) * ydthf_ralsdcp) * zdqsmixdt[jl - 1] + zevaplimmix[jl - 1] = max((zqsmix[jk - 1, jl - 1] - zqx[ncldqv - 1, jk - 1, jl - 1]) / zcorqsmix[jl - 1], + 0.0) + zevaplimliq[jl - 1] = max((zqsliq[jk - 1, jl - 1] - zqx[ncldqv - 1, jk - 1, jl - 1]) / zcorqsliq[jl - 1], + 0.0) + zevaplimice[jl - 1] = max((zqsice[jk - 1, jl - 1] - zqx[ncldqv - 1, jk - 1, jl - 1]) / zcorqsice[jl - 1], + 0.0) + ztmpa = 1.0 / max(za[jk - 1, jl - 1], zepsec) + zliqcld[jl - 1] = zqx[ncldql - 1, jk - 1, jl - 1] * ztmpa + zicecld[jl - 1] = zqx[ncldqi - 1, jk - 1, jl - 1] * ztmpa + zlicld[jl - 1] = zliqcld[jl - 1] + zicecld[jl - 1] + for jl in range(kidia, kfdia + 1): + if zqx[ncldql - 1, jk - 1, jl - 1] < yrecldp_rlmin: + zsolqa[ncldql - 1, ncldqv - 1, jl - 1] = zqx[ncldql - 1, jk - 1, jl - 1] + zsolqa[ncldqv - 1, ncldql - 1, jl - 1] = -zqx[ncldql - 1, jk - 1, jl - 1] + if zqx[ncldqi - 1, jk - 1, jl - 1] < yrecldp_rlmin: + zsolqa[ncldqi - 1, ncldqv - 1, jl - 1] = zqx[ncldqi - 1, jk - 1, jl - 1] + zsolqa[ncldqv - 1, ncldqi - 1, jl - 1] = -zqx[ncldqi - 1, jk - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + zfokoop[jl - 1] = min( + ydthf_rkoop1 - ydthf_rkoop2 * ztp1[jk - 1, jl - 1], + ydthf_r2es * np.exp(ydthf_r3les * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4les)) / + (ydthf_r2es * np.exp(ydthf_r3ies * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4ies)))) + for jl in range(kidia, kfdia + 1): + if ztp1[jk - 1, jl - 1] >= ydcst_rtt or yrecldp_nssopt == 0: + zfac = 1.0 + zfaci = 1.0 + else: + zfac = za[jk - 1, jl - 1] + zfokoop[jl - 1] * (1.0 - za[jk - 1, jl - 1]) + zfaci = ptsphy / yrecldp_rkooptau + if za[jk - 1, jl - 1] > 1.0 - yrecldp_ramin: + zsupsat[jl - 1] = max( + (zqx[ncldqv - 1, jk - 1, jl - 1] - zfac * zqsice[jk - 1, jl - 1]) / zcorqsice[jl - 1], 0.0) + else: + zqp1env = (zqx[ncldqv - 1, jk - 1, jl - 1] - za[jk - 1, jl - 1] * zqsice[jk - 1, jl - 1]) / max( + 1.0 - za[jk - 1, jl - 1], zepsilon) + zsupsat[jl - 1] = max( + (1.0 - za[jk - 1, jl - 1]) * (zqp1env - zfac * zqsice[jk - 1, jl - 1]) / zcorqsice[jl - 1], 0.0) + if zsupsat[jl - 1] > zepsec: + if ztp1[jk - 1, jl - 1] > yrecldp_rthomo: + zsolqa[ncldqv - 1, ncldql - 1, jl - 1] = zsolqa[ncldqv - 1, ncldql - 1, jl - 1] + zsupsat[jl - 1] + zsolqa[ncldql - 1, ncldqv - 1, jl - 1] = zsolqa[ncldql - 1, ncldqv - 1, jl - 1] - zsupsat[jl - 1] + zqxfg[ncldql - 1, jl - 1] = zqxfg[ncldql - 1, jl - 1] + zsupsat[jl - 1] + else: + zsolqa[ncldqv - 1, ncldqi - 1, jl - 1] = zsolqa[ncldqv - 1, ncldqi - 1, jl - 1] + zsupsat[jl - 1] + zsolqa[ncldqi - 1, ncldqv - 1, jl - 1] = zsolqa[ncldqi - 1, ncldqv - 1, jl - 1] - zsupsat[jl - 1] + zqxfg[ncldqi - 1, jl - 1] = zqxfg[ncldqi - 1, jl - 1] + zsupsat[jl - 1] + zsolac[jl - 1] = (1.0 - za[jk - 1, jl - 1]) * zfaci + if psupsat[jk - 1, jl - 1] > zepsec: + if ztp1[jk - 1, jl - 1] > yrecldp_rthomo: + zsolqa[ncldql - 1, ncldql - 1, + jl - 1] = zsolqa[ncldql - 1, ncldql - 1, jl - 1] + psupsat[jk - 1, jl - 1] + zpsupsatsrce[ncldql - 1, jl - 1] = psupsat[jk - 1, jl - 1] + zqxfg[ncldql - 1, jl - 1] = zqxfg[ncldql - 1, jl - 1] + psupsat[jk - 1, jl - 1] + else: + zsolqa[ncldqi - 1, ncldqi - 1, + jl - 1] = zsolqa[ncldqi - 1, ncldqi - 1, jl - 1] + psupsat[jk - 1, jl - 1] + zpsupsatsrce[ncldqi - 1, jl - 1] = psupsat[jk - 1, jl - 1] + zqxfg[ncldqi - 1, jl - 1] = zqxfg[ncldqi - 1, jl - 1] + psupsat[jk - 1, jl - 1] + zsolac[jl - 1] = (1.0 - za[jk - 1, jl - 1]) * zfaci + if jk < klev and jk >= yrecldp_ncldtop: + for jl in range(kidia, kfdia + 1): + plude[jk - 1, jl - 1] = plude[jk - 1, jl - 1] * zdtgdp[jl - 1] + if ldcum[jl - 1] and plude[jk - 1, jl - 1] > yrecldp_rlmin and (plu[jk + 1 - 1, jl - 1] > zepsec): + zsolac[jl - 1] = zsolac[jl - 1] + plude[jk - 1, jl - 1] / plu[jk + 1 - 1, jl - 1] + zalfaw = zfoealfa[jk - 1, jl - 1] + zconvsrce[ncldql - 1, jl - 1] = zalfaw * plude[jk - 1, jl - 1] + zconvsrce[ncldqi - 1, jl - 1] = (1.0 - zalfaw) * plude[jk - 1, jl - 1] + zsolqa[ncldql - 1, ncldql - 1, + jl - 1] = zsolqa[ncldql - 1, ncldql - 1, jl - 1] + zconvsrce[ncldql - 1, jl - 1] + zsolqa[ncldqi - 1, ncldqi - 1, + jl - 1] = zsolqa[ncldqi - 1, ncldqi - 1, jl - 1] + zconvsrce[ncldqi - 1, jl - 1] + else: + plude[jk - 1, jl - 1] = 0.0 + if ldcum[jl - 1]: + zsolqa[ncldqs - 1, ncldqs - 1, + jl - 1] = zsolqa[ncldqs - 1, ncldqs - 1, jl - 1] + psnde[jk - 1, jl - 1] * zdtgdp[jl - 1] + if jk > yrecldp_ncldtop: + for jl in range(kidia, kfdia + 1): + zmf[jl - 1] = max(0.0, (pmfu[jk - 1, jl - 1] + pmfd[jk - 1, jl - 1]) * zdtgdp[jl - 1]) + zacust[jl - 1] = zmf[jl - 1] * zanewm1[jl - 1] + for jm in range(1, nclv + 1): + if not llfall[jm - 1] and iphase[jm - 1] > 0: + for jl in range(kidia, kfdia + 1): + zlcust[jm - 1, jl - 1] = zmf[jl - 1] * zqxnm1[jm - 1, jl - 1] + zconvsrce[jm - 1, jl - 1] = zconvsrce[jm - 1, jl - 1] + zlcust[jm - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + zdtdp = zrdcp * 0.5 * (ztp1[jk - 1 - 1, jl - 1] + ztp1[jk - 1, jl - 1]) / paph[jk - 1, jl - 1] + zdtforc = zdtdp * (pap[jk - 1, jl - 1] - pap[jk - 1 - 1, jl - 1]) + zdqs[jl - 1] = zanewm1[jl - 1] * zdtforc * zdqsmixdt[jl - 1] + for jm in range(1, nclv + 1): + if not llfall[jm - 1] and iphase[jm - 1] > 0: + for jl in range(kidia, kfdia + 1): + zlfinal = max(0.0, zlcust[jm - 1, jl - 1] - zdqs[jl - 1]) + zevap = min(zlcust[jm - 1, jl - 1] - zlfinal, zevaplimmix[jl - 1]) + zlfinal = zlcust[jm - 1, jl - 1] - zevap + zlfinalsum[jl - 1] = zlfinalsum[jl - 1] + zlfinal + zsolqa[jm - 1, jm - 1, jl - 1] = zsolqa[jm - 1, jm - 1, jl - 1] + zlcust[jm - 1, jl - 1] + zsolqa[jm - 1, ncldqv - 1, jl - 1] = zsolqa[jm - 1, ncldqv - 1, jl - 1] + zevap + zsolqa[ncldqv - 1, jm - 1, jl - 1] = zsolqa[ncldqv - 1, jm - 1, jl - 1] - zevap + for jl in range(kidia, kfdia + 1): + if zlfinalsum[jl - 1] < zepsec: + zacust[jl - 1] = 0.0 + zsolac[jl - 1] = zsolac[jl - 1] + zacust[jl - 1] + for jl in range(kidia, kfdia + 1): + if jk < klev: + zmfdn = max(0.0, (pmfu[jk + 1 - 1, jl - 1] + pmfd[jk + 1 - 1, jl - 1]) * zdtgdp[jl - 1]) + zsolab[jl - 1] = zsolab[jl - 1] + zmfdn + zsolqb[ncldql - 1, ncldql - 1, jl - 1] = zsolqb[ncldql - 1, ncldql - 1, jl - 1] + zmfdn + zsolqb[ncldqi - 1, ncldqi - 1, jl - 1] = zsolqb[ncldqi - 1, ncldqi - 1, jl - 1] + zmfdn + zconvsink[ncldql - 1, jl - 1] = zmfdn + zconvsink[ncldqi - 1, jl - 1] = zmfdn + for jl in range(kidia, kfdia + 1): + zldifdt[jl - 1] = yrecldp_rcldiff * ptsphy + if ktype[jl - 1] > 0 and plude[jk - 1, jl - 1] > zepsec: + zldifdt[jl - 1] = yrecldp_rcldiff_convi * zldifdt[jl - 1] + for jl in range(kidia, kfdia + 1): + if zli[jk - 1, jl - 1] > zepsec: + ze = zldifdt[jl - 1] * max(zqsmix[jk - 1, jl - 1] - zqx[ncldqv - 1, jk - 1, jl - 1], 0.0) + zleros = za[jk - 1, jl - 1] * ze + zleros = min(zleros, zevaplimmix[jl - 1]) + zleros = min(zleros, zli[jk - 1, jl - 1]) + zaeros = zleros / zlicld[jl - 1] + zsolac[jl - 1] = zsolac[jl - 1] - zaeros + zsolqa[ncldql - 1, ncldqv - 1, + jl - 1] = zsolqa[ncldql - 1, ncldqv - 1, jl - 1] + zliqfrac[jk - 1, jl - 1] * zleros + zsolqa[ncldqv - 1, ncldql - 1, + jl - 1] = zsolqa[ncldqv - 1, ncldql - 1, jl - 1] - zliqfrac[jk - 1, jl - 1] * zleros + zsolqa[ncldqi - 1, ncldqv - 1, + jl - 1] = zsolqa[ncldqi - 1, ncldqv - 1, jl - 1] + zicefrac[jk - 1, jl - 1] * zleros + zsolqa[ncldqv - 1, ncldqi - 1, + jl - 1] = zsolqa[ncldqv - 1, ncldqi - 1, jl - 1] - zicefrac[jk - 1, jl - 1] * zleros + for jl in range(kidia, kfdia + 1): + zdtdp = zrdcp * ztp1[jk - 1, jl - 1] / pap[jk - 1, jl - 1] + zdpmxdt = zdp[jl - 1] * zqtmst + zmfdn = 0.0 + if jk < klev: + zmfdn = pmfu[jk + 1 - 1, jl - 1] + pmfd[jk + 1 - 1, jl - 1] + zwtot = pvervel[jk - 1, jl - 1] + 0.5 * ydcst_rg * (pmfu[jk - 1, jl - 1] + pmfd[jk - 1, jl - 1] + zmfdn) + zwtot = min(zdpmxdt, max(-zdpmxdt, zwtot)) + zzzdt = phrsw[jk - 1, jl - 1] + phrlw[jk - 1, jl - 1] + zdtdiab = min(zdpmxdt * zdtdp, max(-zdpmxdt * zdtdp, zzzdt)) * ptsphy + ydthf_ralfdcp * zldefr[jl - 1] + zdtforc = zdtdp * zwtot * ptsphy + zdtdiab + zqold[jl - 1] = zqsmix[jk - 1, jl - 1] + ztold[jl - 1] = ztp1[jk - 1, jl - 1] + ztp1[jk - 1, jl - 1] = ztp1[jk - 1, jl - 1] + zdtforc + ztp1[jk - 1, jl - 1] = max(ztp1[jk - 1, jl - 1], 160.0) + llflag[jl - 1] = True + for jl in range(kidia, kfdia + 1): + zqp = 1.0 / pap[jk - 1, jl - 1] + zqsat = ydthf_r2es * (min(1.0, ( + (max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * ydthf_rtwat_rtice_r)**2) * + np.exp(ydthf_r3les * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4les)) + + (1.0 - min(1.0, ( + (max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2)) * np.exp(ydthf_r3ies * + (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4ies))) * zqp + zqsat = min(0.5, zqsat) + zcor = 1.0 / (1.0 - ydcst_retv * zqsat) + zqsat = zqsat * zcor + zcond = (zqsmix[jk - 1, jl - 1] - + zqsat) / (1.0 + zqsat * zcor * + (min(1.0, ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2) * ydthf_r5alvcp * + (1.0 / (ztp1[jk - 1, jl - 1] - ydthf_r4les)**2) + + (1.0 - min(1.0, + ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2)) * ydthf_r5alscp * + (1.0 / (ztp1[jk - 1, jl - 1] - ydthf_r4ies)**2))) + ztp1[jk - 1, jl - 1] = ztp1[jk - 1, jl - 1] + (min(1.0, ( + (max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2) * ydthf_ralvdcp + (1.0 - min(1.0, ( + (max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * ydthf_rtwat_rtice_r)**2)) + * ydthf_ralsdcp) * zcond + zqsmix[jk - 1, jl - 1] = zqsmix[jk - 1, jl - 1] - zcond + zqsat = ydthf_r2es * (min(1.0, ( + (max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * ydthf_rtwat_rtice_r)**2) * + np.exp(ydthf_r3les * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4les)) + + (1.0 - min(1.0, ( + (max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2)) * np.exp(ydthf_r3ies * + (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4ies))) * zqp + zqsat = min(0.5, zqsat) + zcor = 1.0 / (1.0 - ydcst_retv * zqsat) + zqsat = zqsat * zcor + zcond1 = (zqsmix[jk - 1, jl - 1] - + zqsat) / (1.0 + zqsat * zcor * + (min(1.0, ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2) * ydthf_r5alvcp * + (1.0 / (ztp1[jk - 1, jl - 1] - ydthf_r4les)**2) + + (1.0 - min(1.0, + ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2)) * ydthf_r5alscp * + (1.0 / (ztp1[jk - 1, jl - 1] - ydthf_r4ies)**2))) + ztp1[jk - 1, jl - 1] = ztp1[jk - 1, jl - 1] + (min(1.0, ( + (max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2) * ydthf_ralvdcp + (1.0 - min(1.0, ( + (max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * ydthf_rtwat_rtice_r)**2)) + * ydthf_ralsdcp) * zcond1 + zqsmix[jk - 1, jl - 1] = zqsmix[jk - 1, jl - 1] - zcond1 + for jl in range(kidia, kfdia + 1): + zdqs[jl - 1] = zqsmix[jk - 1, jl - 1] - zqold[jl - 1] + zqsmix[jk - 1, jl - 1] = zqold[jl - 1] + ztp1[jk - 1, jl - 1] = ztold[jl - 1] + for jl in range(kidia, kfdia + 1): + if zdqs[jl - 1] > 0.0: + zlevap = za[jk - 1, jl - 1] * min(zdqs[jl - 1], zlicld[jl - 1]) + zlevap = min(zlevap, zevaplimmix[jl - 1]) + zlevap = min(zlevap, max(zqsmix[jk - 1, jl - 1] - zqx[ncldqv - 1, jk - 1, jl - 1], 0.0)) + zlevapl[jl - 1] = zliqfrac[jk - 1, jl - 1] * zlevap + zlevapi[jl - 1] = zicefrac[jk - 1, jl - 1] * zlevap + zsolqa[ncldql - 1, ncldqv - 1, + jl - 1] = zsolqa[ncldql - 1, ncldqv - 1, jl - 1] + zliqfrac[jk - 1, jl - 1] * zlevap + zsolqa[ncldqv - 1, ncldql - 1, + jl - 1] = zsolqa[ncldqv - 1, ncldql - 1, jl - 1] - zliqfrac[jk - 1, jl - 1] * zlevap + zsolqa[ncldqi - 1, ncldqv - 1, + jl - 1] = zsolqa[ncldqi - 1, ncldqv - 1, jl - 1] + zicefrac[jk - 1, jl - 1] * zlevap + zsolqa[ncldqv - 1, ncldqi - 1, + jl - 1] = zsolqa[ncldqv - 1, ncldqi - 1, jl - 1] - zicefrac[jk - 1, jl - 1] * zlevap + for jl in range(kidia, kfdia + 1): + if za[jk - 1, jl - 1] > zepsec and zdqs[jl - 1] <= -yrecldp_rlmin: + zlcond1[jl - 1] = max(-zdqs[jl - 1], 0.0) + if za[jk - 1, jl - 1] > 0.99: + zcor = 1.0 / (1.0 - ydcst_retv * zqsmix[jk - 1, jl - 1]) + zcdmax = (zqx[ncldqv - 1, jk - 1, jl - 1] - zqsmix[jk - 1, jl - 1]) / ( + 1.0 + zcor * zqsmix[jk - 1, jl - 1] * + (min(1.0, ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2) * ydthf_r5alvcp * + (1.0 / (ztp1[jk - 1, jl - 1] - ydthf_r4les)**2) + + (1.0 - min(1.0, ((max(ydthf_rtice, min(ydthf_rtwat, ztp1[jk - 1, jl - 1])) - ydthf_rtice) * + ydthf_rtwat_rtice_r)**2)) * ydthf_r5alscp * + (1.0 / (ztp1[jk - 1, jl - 1] - ydthf_r4ies)**2))) + else: + zcdmax = (zqx[ncldqv - 1, jk - 1, jl - 1] - + za[jk - 1, jl - 1] * zqsmix[jk - 1, jl - 1]) / za[jk - 1, jl - 1] + zlcond1[jl - 1] = max(min(zlcond1[jl - 1], zcdmax), 0.0) + zlcond1[jl - 1] = za[jk - 1, jl - 1] * zlcond1[jl - 1] + if zlcond1[jl - 1] < yrecldp_rlmin: + zlcond1[jl - 1] = 0.0 + if ztp1[jk - 1, jl - 1] > yrecldp_rthomo: + zsolqa[ncldqv - 1, ncldql - 1, jl - 1] = zsolqa[ncldqv - 1, ncldql - 1, jl - 1] + zlcond1[jl - 1] + zsolqa[ncldql - 1, ncldqv - 1, jl - 1] = zsolqa[ncldql - 1, ncldqv - 1, jl - 1] - zlcond1[jl - 1] + zqxfg[ncldql - 1, jl - 1] = zqxfg[ncldql - 1, jl - 1] + zlcond1[jl - 1] + else: + zsolqa[ncldqv - 1, ncldqi - 1, jl - 1] = zsolqa[ncldqv - 1, ncldqi - 1, jl - 1] + zlcond1[jl - 1] + zsolqa[ncldqi - 1, ncldqv - 1, jl - 1] = zsolqa[ncldqi - 1, ncldqv - 1, jl - 1] - zlcond1[jl - 1] + zqxfg[ncldqi - 1, jl - 1] = zqxfg[ncldqi - 1, jl - 1] + zlcond1[jl - 1] + for jl in range(kidia, kfdia + 1): + if zdqs[jl - 1] <= -yrecldp_rlmin and za[jk - 1, jl - 1] < 1.0 - zepsec: + zsigk = pap[jk - 1, jl - 1] / paph[klev + 1 - 1, jl - 1] + if zsigk > 0.8: + zrhc = yrecldp_ramid + (1.0 - yrecldp_ramid) * ((zsigk - 0.8) / 0.2)**2 + else: + zrhc = yrecldp_ramid + if yrecldp_nssopt == 0: + zqe = (zqx[ncldqv - 1, jk - 1, jl - 1] - za[jk - 1, jl - 1] * zqsice[jk - 1, jl - 1]) / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zqe = max(0.0, zqe) + elif yrecldp_nssopt == 1: + zqe = (zqx[ncldqv - 1, jk - 1, jl - 1] - za[jk - 1, jl - 1] * zqsice[jk - 1, jl - 1]) / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zqe = max(0.0, zqe) + elif yrecldp_nssopt == 2: + zqe = zqx[ncldqv - 1, jk - 1, jl - 1] + elif yrecldp_nssopt == 3: + zqe = zqx[ncldqv - 1, jk - 1, jl - 1] + zli[jk - 1, jl - 1] + if ztp1[jk - 1, jl - 1] >= ydcst_rtt or yrecldp_nssopt == 0: + zfac = 1.0 + else: + zfac = zfokoop[jl - 1] + if zqe >= zrhc * zqsice[jk - 1, jl - 1] * zfac and zqe < zqsice[jk - 1, jl - 1] * zfac: + zacond = -(1.0 - za[jk - 1, jl - 1]) * zfac * zdqs[jl - 1] / max( + 2.0 * (zfac * zqsice[jk - 1, jl - 1] - zqe), zepsec) + zacond = min(zacond, 1.0 - za[jk - 1, jl - 1]) + zlcond2[jl - 1] = -zfac * zdqs[jl - 1] * 0.5 * zacond + zzdl = 2.0 * (zfac * zqsice[jk - 1, jl - 1] - zqe) / max(zepsec, 1.0 - za[jk - 1, jl - 1]) + if zfac * zdqs[jl - 1] < -zzdl: + zlcondlim = (za[jk - 1, jl - 1] - + 1.0) * zfac * zdqs[jl - 1] - zfac * zqsice[jk - 1, jl - 1] + zqx[ncldqv - 1, + jk - 1, jl - 1] + zlcond2[jl - 1] = min(zlcond2[jl - 1], zlcondlim) + zlcond2[jl - 1] = max(zlcond2[jl - 1], 0.0) + if zlcond2[jl - 1] < yrecldp_rlmin or 1.0 - za[jk - 1, jl - 1] < zepsec: + zlcond2[jl - 1] = 0.0 + zacond = 0.0 + if zlcond2[jl - 1] == 0.0: + zacond = 0.0 + zsolac[jl - 1] = zsolac[jl - 1] + zacond + if ztp1[jk - 1, jl - 1] > yrecldp_rthomo: + zsolqa[ncldqv - 1, ncldql - 1, + jl - 1] = zsolqa[ncldqv - 1, ncldql - 1, jl - 1] + zlcond2[jl - 1] + zsolqa[ncldql - 1, ncldqv - 1, + jl - 1] = zsolqa[ncldql - 1, ncldqv - 1, jl - 1] - zlcond2[jl - 1] + zqxfg[ncldql - 1, jl - 1] = zqxfg[ncldql - 1, jl - 1] + zlcond2[jl - 1] + else: + zsolqa[ncldqv - 1, ncldqi - 1, + jl - 1] = zsolqa[ncldqv - 1, ncldqi - 1, jl - 1] + zlcond2[jl - 1] + zsolqa[ncldqi - 1, ncldqv - 1, + jl - 1] = zsolqa[ncldqi - 1, ncldqv - 1, jl - 1] - zlcond2[jl - 1] + zqxfg[ncldqi - 1, jl - 1] = zqxfg[ncldqi - 1, jl - 1] + zlcond2[jl - 1] + if idepice == 1: + for jl in range(kidia, kfdia + 1): + if za[jk - 1 - 1, jl - 1] < yrecldp_rcldtopcf and za[jk - 1, jl - 1] >= yrecldp_rcldtopcf: + zcldtopdist[jl - 1] = 0.0 + else: + zcldtopdist[jl - 1] = zcldtopdist[jl - 1] + zdp[jl - 1] / (zrho[jl - 1] * ydcst_rg) + if ztp1[jk - 1, jl - 1] < ydcst_rtt and zqxfg[ncldql - 1, jl - 1] > yrecldp_rlmin: + zvpice = ydthf_r2es * np.exp(ydthf_r3ies * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4ies)) * ydcst_rv / ydcst_rd + zvpliq = zvpice * zfokoop[jl - 1] + zicenuclei[jl - 1] = 1000.0 * np.exp(12.96 * (zvpliq - zvpice) / zvpliq - 0.639) + zadd = ydcst_rlstt * (ydcst_rlstt / + (ydcst_rv * ztp1[jk - 1, jl - 1]) - 1.0) / (0.024 * ztp1[jk - 1, jl - 1]) + zbdd = ydcst_rv * ztp1[jk - 1, jl - 1] * pap[jk - 1, jl - 1] / (2.21 * zvpice) + zcvds = 7.8 * (zicenuclei[jl - 1] / + zrho[jl - 1])**0.666 * (zvpliq - zvpice) / (8.87 * (zadd + zbdd) * zvpice) + zice0 = max(zicecld[jl - 1], zicenuclei[jl - 1] * yrecldp_riceinit / zrho[jl - 1]) + zinew = (0.666 * zcvds * ptsphy + zice0**0.666)**1.5 + zdepos = max(za[jk - 1, jl - 1] * (zinew - zice0), 0.0) + zdepos = min(zdepos, zqxfg[ncldql - 1, jl - 1]) + zinfactor = min(zicenuclei[jl - 1] / 15000.0, 1.0) + zdepos = zdepos * min( + zinfactor + (1.0 - zinfactor) * + (yrecldp_rdepliqrefrate + zcldtopdist[jl - 1] / yrecldp_rdepliqrefdepth), 1.0) + zsolqa[ncldql - 1, ncldqi - 1, jl - 1] = zsolqa[ncldql - 1, ncldqi - 1, jl - 1] + zdepos + zsolqa[ncldqi - 1, ncldql - 1, jl - 1] = zsolqa[ncldqi - 1, ncldql - 1, jl - 1] - zdepos + zqxfg[ncldqi - 1, jl - 1] = zqxfg[ncldqi - 1, jl - 1] + zdepos + zqxfg[ncldql - 1, jl - 1] = zqxfg[ncldql - 1, jl - 1] - zdepos + elif idepice == 2: + for jl in range(kidia, kfdia + 1): + if za[jk - 1 - 1, jl - 1] < yrecldp_rcldtopcf and za[jk - 1, jl - 1] >= yrecldp_rcldtopcf: + zcldtopdist[jl - 1] = 0.0 + else: + zcldtopdist[jl - 1] = zcldtopdist[jl - 1] + zdp[jl - 1] / (zrho[jl - 1] * ydcst_rg) + if ztp1[jk - 1, jl - 1] < ydcst_rtt and zqxfg[ncldql - 1, jl - 1] > yrecldp_rlmin: + zvpice = ydthf_r2es * np.exp(ydthf_r3ies * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4ies)) * ydcst_rv / ydcst_rd + zvpliq = zvpice * zfokoop[jl - 1] + zicenuclei[jl - 1] = 1000.0 * np.exp(12.96 * (zvpliq - zvpice) / zvpliq - 0.639) + zice0 = max(zicecld[jl - 1], zicenuclei[jl - 1] * yrecldp_riceinit / zrho[jl - 1]) + ztcg = 1.0 + zfacx1i = 1.0 + zaplusb = yrecldp_rcl_apb1 * zvpice - yrecldp_rcl_apb2 * zvpice * ztp1[jk - 1, jl - 1] + pap[ + jk - 1, jl - 1] * yrecldp_rcl_apb3 * ztp1[jk - 1, jl - 1]**3.0 + zcorrfac = (1.0 / zrho[jl - 1])**0.5 + zcorrfac2 = (ztp1[jk - 1, jl - 1] / 273.0)**1.5 * (393.0 / (ztp1[jk - 1, jl - 1] + 120.0)) + zpr02 = zrho[jl - 1] * zice0 * yrecldp_rcl_const1i / (ztcg * zfacx1i) + zterm1 = (zvpliq - + zvpice) * ztp1[jk - 1, jl - + 1]**2.0 * zvpice * zcorrfac2 * ztcg * yrecldp_rcl_const2i * zfacx1i / ( + zrho[jl - 1] * zaplusb * zvpice) + zterm2 = 0.65 * yrecldp_rcl_const6i * zpr02**yrecldp_rcl_const4i + yrecldp_rcl_const3i * zcorrfac**0.5 * zrho[ + jl - 1]**0.5 * zpr02**yrecldp_rcl_const5i / zcorrfac2**0.5 + zdepos = max(za[jk - 1, jl - 1] * zterm1 * zterm2 * ptsphy, 0.0) + zdepos = min(zdepos, zqxfg[ncldql - 1, jl - 1]) + zinfactor = min(zicenuclei[jl - 1] / 15000.0, 1.0) + zdepos = zdepos * min( + zinfactor + (1.0 - zinfactor) * + (yrecldp_rdepliqrefrate + zcldtopdist[jl - 1] / yrecldp_rdepliqrefdepth), 1.0) + zsolqa[ncldql - 1, ncldqi - 1, jl - 1] = zsolqa[ncldql - 1, ncldqi - 1, jl - 1] + zdepos + zsolqa[ncldqi - 1, ncldql - 1, jl - 1] = zsolqa[ncldqi - 1, ncldql - 1, jl - 1] - zdepos + zqxfg[ncldqi - 1, jl - 1] = zqxfg[ncldqi - 1, jl - 1] + zdepos + zqxfg[ncldql - 1, jl - 1] = zqxfg[ncldql - 1, jl - 1] - zdepos + for jl in range(kidia, kfdia + 1): + ztmpa = 1.0 / max(za[jk - 1, jl - 1], zepsec) + zliqcld[jl - 1] = zqxfg[ncldql - 1, jl - 1] * ztmpa + zicecld[jl - 1] = zqxfg[ncldqi - 1, jl - 1] * ztmpa + zlicld[jl - 1] = zliqcld[jl - 1] + zicecld[jl - 1] + for jm in range(1, nclv + 1): + if llfall[jm - 1] or jm == ncldqi: + for jl in range(kidia, kfdia + 1): + if jk > yrecldp_ncldtop: + zfallsrce[jm - 1, jl - 1] = zpfplsx[jm - 1, jk - 1, jl - 1] * zdtgdp[jl - 1] + zsolqa[jm - 1, jm - 1, jl - 1] = zsolqa[jm - 1, jm - 1, jl - 1] + zfallsrce[jm - 1, jl - 1] + zqxfg[jm - 1, jl - 1] = zqxfg[jm - 1, jl - 1] + zfallsrce[jm - 1, jl - 1] + zqpretot[jl - 1] = zqpretot[jl - 1] + zqxfg[jm - 1, jl - 1] + if yrecldp_laericesed and jm == ncldqi: + zre_ice = pre_ice[jk - 1, jl - 1] + zvqx[ncldqi - 1] = 0.002 * zre_ice**1.0 + zfall = zvqx[jm - 1] * zrho[jl - 1] + zfallsink[jm - 1, jl - 1] = zdtgdp[jl - 1] * zfall + for jl in range(kidia, kfdia + 1): + if zqpretot[jl - 1] > zepsec: + zcovptot[jl - 1] = 1.0 - (1.0 - zcovptot[jl - 1]) * (1.0 - max( + za[jk - 1, jl - 1], za[jk - 1 - 1, jl - 1])) / (1.0 - min(za[jk - 1 - 1, jl - 1], 1.0 - 1e-06)) + zcovptot[jl - 1] = max(zcovptot[jl - 1], yrecldp_rcovpmin) + zcovpclr[jl - 1] = max(0.0, zcovptot[jl - 1] - za[jk - 1, jl - 1]) + zraincld[jl - 1] = zqxfg[ncldqr - 1, jl - 1] / zcovptot[jl - 1] + zsnowcld[jl - 1] = zqxfg[ncldqs - 1, jl - 1] / zcovptot[jl - 1] + zcovpmax[jl - 1] = max(zcovptot[jl - 1], zcovpmax[jl - 1]) + else: + zraincld[jl - 1] = 0.0 + zsnowcld[jl - 1] = 0.0 + zcovptot[jl - 1] = 0.0 + zcovpclr[jl - 1] = 0.0 + zcovpmax[jl - 1] = 0.0 + for jl in range(kidia, kfdia + 1): + if ztp1[jk - 1, jl - 1] <= ydcst_rtt: + if zicecld[jl - 1] > zepsec: + zzco = ptsphy * yrecldp_rsnowlin1 * np.exp(yrecldp_rsnowlin2 * (ztp1[jk - 1, jl - 1] - ydcst_rtt)) + if yrecldp_laericeauto: + zlcrit = picrit_aer[jk - 1, jl - 1] + zzco = zzco * (yrecldp_rnice / pnice[jk - 1, jl - 1])**0.333 + else: + zlcrit = yrecldp_rlcritsnow + zsnowaut[jl - 1] = zzco * (1.0 - np.exp(-(zicecld[jl - 1] / zlcrit)**2)) + zsolqb[ncldqi - 1, ncldqs - 1, jl - 1] = zsolqb[ncldqi - 1, ncldqs - 1, jl - 1] + zsnowaut[jl - 1] + if zliqcld[jl - 1] > zepsec: + if iwarmrain == 1: + zzco = yrecldp_rkconv * ptsphy + if yrecldp_laerliqautolsp: + zlcrit = plcrit_aer[jk - 1, jl - 1] + zzco = zzco * (yrecldp_rccn / pccn[jk - 1, jl - 1])**0.333 + elif plsm[jl - 1] > 0.5: + zlcrit = yrecldp_rclcrit_land + else: + zlcrit = yrecldp_rclcrit_sea + zprecip = (zpfplsx[ncldqs - 1, jk - 1, jl - 1] + zpfplsx[ncldqr - 1, jk - 1, jl - 1]) / max( + zepsec, zcovptot[jl - 1]) + zcfpr = 1.0 + yrecldp_rprc1 * np.sqrt(max(zprecip, 0.0)) + if yrecldp_laerliqcoll: + zcfpr = zcfpr * (yrecldp_rccn / pccn[jk - 1, jl - 1])**0.333 + zzco = zzco * zcfpr + zlcrit = zlcrit / max(zcfpr, zepsec) + if zliqcld[jl - 1] / zlcrit < 20.0: + zrainaut[jl - 1] = zzco * (1.0 - np.exp(-(zliqcld[jl - 1] / zlcrit)**2)) + else: + zrainaut[jl - 1] = zzco + if ztp1[jk - 1, jl - 1] <= ydcst_rtt: + zsolqb[ncldql - 1, ncldqs - 1, + jl - 1] = zsolqb[ncldql - 1, ncldqs - 1, jl - 1] + zrainaut[jl - 1] + else: + zsolqb[ncldql - 1, ncldqr - 1, + jl - 1] = zsolqb[ncldql - 1, ncldqr - 1, jl - 1] + zrainaut[jl - 1] + elif iwarmrain == 2: + if plsm[jl - 1] > 0.5: + zconst = yrecldp_rcl_kk_cloud_num_land + zlcrit = yrecldp_rclcrit_land + else: + zconst = yrecldp_rcl_kk_cloud_num_sea + zlcrit = yrecldp_rclcrit_sea + if zliqcld[jl - 1] > zlcrit: + zrainaut[jl - 1] = 1.5 * za[jk - 1, jl - 1] * ptsphy * yrecldp_rcl_kkaau * zliqcld[ + jl - 1]**yrecldp_rcl_kkbauq * zconst**yrecldp_rcl_kkbaun + zrainaut[jl - 1] = min(zrainaut[jl - 1], zqxfg[ncldql - 1, jl - 1]) + if zrainaut[jl - 1] < zepsec: + zrainaut[jl - 1] = 0.0 + zrainacc[jl - 1] = 2.0 * za[jk - 1, jl - 1] * ptsphy * yrecldp_rcl_kkaac * ( + zliqcld[jl - 1] * zraincld[jl - 1])**yrecldp_rcl_kkbac + zrainacc[jl - 1] = min(zrainacc[jl - 1], zqxfg[ncldql - 1, jl - 1]) + if zrainacc[jl - 1] < zepsec: + zrainacc[jl - 1] = 0.0 + else: + zrainaut[jl - 1] = 0.0 + zrainacc[jl - 1] = 0.0 + if ztp1[jk - 1, jl - 1] <= ydcst_rtt: + zsolqa[ncldql - 1, ncldqs - 1, + jl - 1] = zsolqa[ncldql - 1, ncldqs - 1, jl - 1] + zrainaut[jl - 1] + zsolqa[ncldql - 1, ncldqs - 1, + jl - 1] = zsolqa[ncldql - 1, ncldqs - 1, jl - 1] + zrainacc[jl - 1] + zsolqa[ncldqs - 1, ncldql - 1, + jl - 1] = zsolqa[ncldqs - 1, ncldql - 1, jl - 1] - zrainaut[jl - 1] + zsolqa[ncldqs - 1, ncldql - 1, + jl - 1] = zsolqa[ncldqs - 1, ncldql - 1, jl - 1] - zrainacc[jl - 1] + else: + zsolqa[ncldql - 1, ncldqr - 1, + jl - 1] = zsolqa[ncldql - 1, ncldqr - 1, jl - 1] + zrainaut[jl - 1] + zsolqa[ncldql - 1, ncldqr - 1, + jl - 1] = zsolqa[ncldql - 1, ncldqr - 1, jl - 1] + zrainacc[jl - 1] + zsolqa[ncldqr - 1, ncldql - 1, + jl - 1] = zsolqa[ncldqr - 1, ncldql - 1, jl - 1] - zrainaut[jl - 1] + zsolqa[ncldqr - 1, ncldql - 1, + jl - 1] = zsolqa[ncldqr - 1, ncldql - 1, jl - 1] - zrainacc[jl - 1] + if iwarmrain > 1: + for jl in range(kidia, kfdia + 1): + if ztp1[jk - 1, jl - 1] <= ydcst_rtt and zliqcld[jl - 1] > zepsec: + zfallcorr = (yrecldp_rdensref / zrho[jl - 1])**0.4 + if zsnowcld[jl - 1] > zepsec and zcovptot[jl - 1] > 0.01: + zsnowrime[jl - 1] = 0.3 * zcovptot[jl - 1] * ptsphy * yrecldp_rcl_const7s * zfallcorr * ( + zrho[jl - 1] * zsnowcld[jl - 1] * yrecldp_rcl_const1s)**yrecldp_rcl_const8s + zsnowrime[jl - 1] = min(zsnowrime[jl - 1], 1.0) + zsolqb[ncldql - 1, ncldqs - 1, + jl - 1] = zsolqb[ncldql - 1, ncldqs - 1, jl - 1] + zsnowrime[jl - 1] + for jl in range(kidia, kfdia + 1): + zicetot[jl - 1] = zqxfg[ncldqi - 1, jl - 1] + zqxfg[ncldqs - 1, jl - 1] + zmeltmax[jl - 1] = 0.0 + if zicetot[jl - 1] > zepsec and ztp1[jk - 1, jl - 1] > ydcst_rtt: + zsubsat = max(zqsice[jk - 1, jl - 1] - zqx[ncldqv - 1, jk - 1, jl - 1], 0.0) + ztdmtw0 = ztp1[jk - 1, jl - 1] - ydcst_rtt - zsubsat * (ztw1 + ztw2 * + (pap[jk - 1, jl - 1] - ztw3) - ztw4 * + (ztp1[jk - 1, jl - 1] - ztw5)) + zcons1 = abs(ptsphy * (1.0 + 0.5 * ztdmtw0) / yrecldp_rtaumel) + zmeltmax[jl - 1] = max(ztdmtw0 * zcons1 * zrldcp, 0.0) + for jm in range(1, nclv + 1): + if iphase[jm - 1] == 2: + for jl in range(kidia, kfdia + 1): + if zmeltmax[jl - 1] > zepsec and zicetot[jl - 1] > zepsec: + zalfa2 = zqxfg[jm - 1, jl - 1] / zicetot[jl - 1] + zmelt = min(zqxfg[jm - 1, jl - 1], zalfa2 * zmeltmax[jl - 1]) + zqxfg[jm - 1, jl - 1] = zqxfg[jm - 1, jl - 1] - zmelt + zqxfg[imelt[jm - 1] - 1, jl - 1] = zqxfg[imelt[jm - 1] - 1, jl - 1] + zmelt + zsolqa[jm - 1, imelt[jm - 1] - 1, jl - 1] = zsolqa[jm - 1, imelt[jm - 1] - 1, jl - 1] + zmelt + zsolqa[imelt[jm - 1] - 1, jm - 1, jl - 1] = zsolqa[imelt[jm - 1] - 1, jm - 1, jl - 1] - zmelt + for jl in range(kidia, kfdia + 1): + if zqx[ncldqr - 1, jk - 1, jl - 1] > zepsec: + if ztp1[jk - 1, jl - 1] <= ydcst_rtt and ztp1[jk - 1 - 1, jl - 1] > ydcst_rtt: + zqpretot[jl - 1] = max(zqx[ncldqs - 1, jk - 1, jl - 1] + zqx[ncldqr - 1, jk - 1, jl - 1], zepsec) + prainfrac_toprfz[jl - 1] = zqx[ncldqr - 1, jk - 1, jl - 1] / zqpretot[jl - 1] + if prainfrac_toprfz[jl - 1] > 0.8: + llrainliq[jl - 1] = True + else: + llrainliq[jl - 1] = False + if ztp1[jk - 1, jl - 1] < ydcst_rtt: + if prainfrac_toprfz[jl - 1] > 0.8: + zlambda = (yrecldp_rcl_fac1 / + (zrho[jl - 1] * zqx[ncldqr - 1, jk - 1, jl - 1]))**yrecldp_rcl_fac2 + ztemp = yrecldp_rcl_fzrab * (ztp1[jk - 1, jl - 1] - ydcst_rtt) + zfrz = ptsphy * (yrecldp_rcl_const5r / zrho[jl - 1]) * (np.exp(ztemp) - + 1.0) * zlambda**yrecldp_rcl_const6r + zfrzmax[jl - 1] = max(zfrz, 0.0) + else: + zcons1 = abs(ptsphy * (1.0 + 0.5 * (ydcst_rtt - ztp1[jk - 1, jl - 1])) / yrecldp_rtaumel) + zfrzmax[jl - 1] = max((ydcst_rtt - ztp1[jk - 1, jl - 1]) * zcons1 * zrldcp, 0.0) + if zfrzmax[jl - 1] > zepsec: + zfrz = min(zqx[ncldqr - 1, jk - 1, jl - 1], zfrzmax[jl - 1]) + zsolqa[ncldqr - 1, ncldqs - 1, jl - 1] = zsolqa[ncldqr - 1, ncldqs - 1, jl - 1] + zfrz + zsolqa[ncldqs - 1, ncldqr - 1, jl - 1] = zsolqa[ncldqs - 1, ncldqr - 1, jl - 1] - zfrz + for jl in range(kidia, kfdia + 1): + zfrzmax[jl - 1] = max((yrecldp_rthomo - ztp1[jk - 1, jl - 1]) * zrldcp, 0.0) + + for jl in range(kidia, kfdia + 1): + if zfrzmax[jl - 1] > zepsec and zqxfg[ncldql - 1, jl - 1] > zepsec: + zfrz = min(zqxfg[ncldql - 1, jl - 1], zfrzmax[jl - 1]) + zsolqa[ncldql - 1, imelt[ncldql - 1] - 1, + jl - 1] = zsolqa[ncldql - 1, imelt[ncldql - 1] - 1, jl - 1] + zfrz + zsolqa[imelt[ncldql - 1] - 1, ncldql - 1, + jl - 1] = zsolqa[imelt[ncldql - 1] - 1, ncldql - 1, jl - 1] - zfrz + if ievaprain == 1: + for jl in range(kidia, kfdia + 1): + zzrh = yrecldp_rprecrhmax + (1.0 - yrecldp_rprecrhmax) * zcovpmax[jl - 1] / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zzrh = min(max(zzrh, yrecldp_rprecrhmax), 1.0) + zqe = (zqx[ncldqv - 1, jk - 1, jl - 1] - za[jk - 1, jl - 1] * zqsliq[jk - 1, jl - 1]) / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zqe = max(0.0, min(zqe, zqsliq[jk - 1, jl - 1])) + llo1 = zcovpclr[jl - 1] > zepsec and zqxfg[ncldqr - 1, + jl - 1] > zepsec and (zqe < zzrh * zqsliq[jk - 1, jl - 1]) + if llo1: + zpreclr = zqxfg[ncldqr - 1, jl - 1] * zcovpclr[jl - 1] / (max( + abs(zcovptot[jl - 1] * zdtgdp[jl - 1]), zepsilon) * np.sign(zcovptot[jl - 1] * zdtgdp[jl - 1])) + zbeta1 = np.sqrt( + pap[jk - 1, jl - 1] / paph[klev + 1 - 1, jl - 1]) / yrecldp_rvrfactor * zpreclr / max( + zcovpclr[jl - 1], zepsec) + zbeta = ydcst_rg * yrecldp_rpecons * 0.5 * zbeta1**0.5777 + zdenom = 1.0 + zbeta * ptsphy * zcorqsliq[jl - 1] + zdpr = zcovpclr[jl - 1] * zbeta * (zqsliq[jk - 1, jl - 1] - zqe) / zdenom * zdp[jl - 1] * zrg_r + zdpevap = zdpr * zdtgdp[jl - 1] + zevap = min(zdpevap, zqxfg[ncldqr - 1, jl - 1]) + zsolqa[ncldqr - 1, ncldqv - 1, jl - 1] = zsolqa[ncldqr - 1, ncldqv - 1, jl - 1] + zevap + zsolqa[ncldqv - 1, ncldqr - 1, jl - 1] = zsolqa[ncldqv - 1, ncldqr - 1, jl - 1] - zevap + zcovptot[jl - 1] = max( + yrecldp_rcovpmin, zcovptot[jl - 1] - + max(0.0, (zcovptot[jl - 1] - za[jk - 1, jl - 1]) * zevap / zqxfg[ncldqr - 1, jl - 1])) + zqxfg[ncldqr - 1, jl - 1] = zqxfg[ncldqr - 1, jl - 1] - zevap + elif ievaprain == 2: + for jl in range(kidia, kfdia + 1): + zzrh = yrecldp_rprecrhmax + (1.0 - yrecldp_rprecrhmax) * zcovpmax[jl - 1] / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zzrh = min(max(zzrh, yrecldp_rprecrhmax), 1.0) + zzrh = min(0.8, zzrh) + zqe = max(0.0, min(zqx[ncldqv - 1, jk - 1, jl - 1], zqsliq[jk - 1, jl - 1])) + llo1 = zcovpclr[jl - 1] > zepsec and zqxfg[ncldqr - 1, + jl - 1] > zepsec and (zqe < zzrh * zqsliq[jk - 1, jl - 1]) + if llo1: + zpreclr = zqxfg[ncldqr - 1, jl - 1] / zcovptot[jl - 1] + zfallcorr = (yrecldp_rdensref / zrho[jl - 1])**0.4 + zesatliq = ydcst_rv / ydcst_rd * (ydthf_r2es * np.exp(ydthf_r3les * + (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4les))) + zlambda = (yrecldp_rcl_fac1 / (zrho[jl - 1] * zpreclr))**yrecldp_rcl_fac2 + zevap_denom = yrecldp_rcl_cdenom1 * zesatliq - yrecldp_rcl_cdenom2 * ztp1[ + jk - 1, jl - 1] * zesatliq + yrecldp_rcl_cdenom3 * ztp1[jk - 1, jl - 1]**3.0 * pap[jk - 1, + jl - 1] + zcorr2 = (ztp1[jk - 1, jl - 1] / 273.0)**1.5 * 393.0 / (ztp1[jk - 1, jl - 1] + 120.0) + zka = yrecldp_rcl_ka273 * zcorr2 + zsubsat = max(zzrh * zqsliq[jk - 1, jl - 1] - zqe, 0.0) + zbeta = 0.5 / zqsliq[jk - 1, + jl - 1] * ztp1[jk - 1, jl - 1]**2.0 * zesatliq * yrecldp_rcl_const1r * ( + zcorr2 / + zevap_denom) * (0.78 / zlambda**yrecldp_rcl_const4r + yrecldp_rcl_const2r * + (zrho[jl - 1] * zfallcorr)**0.5 / + (zcorr2**0.5 * zlambda**yrecldp_rcl_const3r)) + zdenom = 1.0 + zbeta * ptsphy + zdpevap = zcovpclr[jl - 1] * zbeta * ptsphy * zsubsat / zdenom + zevap = min(zdpevap, zqxfg[ncldqr - 1, jl - 1]) + zsolqa[ncldqr - 1, ncldqv - 1, jl - 1] = zsolqa[ncldqr - 1, ncldqv - 1, jl - 1] + zevap + zsolqa[ncldqv - 1, ncldqr - 1, jl - 1] = zsolqa[ncldqv - 1, ncldqr - 1, jl - 1] - zevap + zcovptot[jl - 1] = max( + yrecldp_rcovpmin, zcovptot[jl - 1] - + max(0.0, (zcovptot[jl - 1] - za[jk - 1, jl - 1]) * zevap / zqxfg[ncldqr - 1, jl - 1])) + zqxfg[ncldqr - 1, jl - 1] = zqxfg[ncldqr - 1, jl - 1] - zevap + if ievapsnow == 1: + for jl in range(kidia, kfdia + 1): + zzrh = yrecldp_rprecrhmax + (1.0 - yrecldp_rprecrhmax) * zcovpmax[jl - 1] / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zzrh = min(max(zzrh, yrecldp_rprecrhmax), 1.0) + zqe = (zqx[ncldqv - 1, jk - 1, jl - 1] - za[jk - 1, jl - 1] * zqsice[jk - 1, jl - 1]) / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zqe = max(0.0, min(zqe, zqsice[jk - 1, jl - 1])) + llo1 = zcovpclr[jl - 1] > zepsec and zqxfg[ncldqs - 1, + jl - 1] > zepsec and (zqe < zzrh * zqsice[jk - 1, jl - 1]) + if llo1: + zpreclr = zqxfg[ncldqs - 1, jl - 1] * zcovpclr[jl - 1] / (max( + abs(zcovptot[jl - 1] * zdtgdp[jl - 1]), zepsilon) * np.sign(zcovptot[jl - 1] * zdtgdp[jl - 1])) + zbeta1 = np.sqrt( + pap[jk - 1, jl - 1] / paph[klev + 1 - 1, jl - 1]) / yrecldp_rvrfactor * zpreclr / max( + zcovpclr[jl - 1], zepsec) + zbeta = ydcst_rg * yrecldp_rpecons * zbeta1**0.5777 + zdenom = 1.0 + zbeta * ptsphy * zcorqsice[jl - 1] + zdpr = zcovpclr[jl - 1] * zbeta * (zqsice[jk - 1, jl - 1] - zqe) / zdenom * zdp[jl - 1] * zrg_r + zdpevap = zdpr * zdtgdp[jl - 1] + zevap = min(zdpevap, zqxfg[ncldqs - 1, jl - 1]) + zsolqa[ncldqs - 1, ncldqv - 1, jl - 1] = zsolqa[ncldqs - 1, ncldqv - 1, jl - 1] + zevap + zsolqa[ncldqv - 1, ncldqs - 1, jl - 1] = zsolqa[ncldqv - 1, ncldqs - 1, jl - 1] - zevap + zcovptot[jl - 1] = max( + yrecldp_rcovpmin, zcovptot[jl - 1] - + max(0.0, (zcovptot[jl - 1] - za[jk - 1, jl - 1]) * zevap / zqxfg[ncldqs - 1, jl - 1])) + zqxfg[ncldqs - 1, jl - 1] = zqxfg[ncldqs - 1, jl - 1] - zevap + elif ievapsnow == 2: + for jl in range(kidia, kfdia + 1): + zzrh = yrecldp_rprecrhmax + (1.0 - yrecldp_rprecrhmax) * zcovpmax[jl - 1] / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zzrh = min(max(zzrh, yrecldp_rprecrhmax), 1.0) + zqe = (zqx[ncldqv - 1, jk - 1, jl - 1] - za[jk - 1, jl - 1] * zqsice[jk - 1, jl - 1]) / max( + zepsec, 1.0 - za[jk - 1, jl - 1]) + zqe = max(0.0, min(zqe, zqsice[jk - 1, jl - 1])) + llo1 = zcovpclr[jl - 1] > zepsec and zqx[ncldqs - 1, jk - 1, + jl - 1] > zepsec and (zqe < zzrh * zqsice[jk - 1, jl - 1]) + if llo1: + zpreclr = zqx[ncldqs - 1, jk - 1, jl - 1] / zcovptot[jl - 1] + zvpice = ydthf_r2es * np.exp(ydthf_r3ies * (ztp1[jk - 1, jl - 1] - ydcst_rtt) / + (ztp1[jk - 1, jl - 1] - ydthf_r4ies)) * ydcst_rv / ydcst_rd + ztcg = 1.0 + zfacx1s = 1.0 + zaplusb = yrecldp_rcl_apb1 * zvpice - yrecldp_rcl_apb2 * zvpice * ztp1[jk - 1, jl - 1] + pap[ + jk - 1, jl - 1] * yrecldp_rcl_apb3 * ztp1[jk - 1, jl - 1]**3 + zcorrfac = (1.0 / zrho[jl - 1])**0.5 + zcorrfac2 = (ztp1[jk - 1, jl - 1] / 273.0)**1.5 * (393.0 / (ztp1[jk - 1, jl - 1] + 120.0)) + zpr02 = zrho[jl - 1] * zpreclr * yrecldp_rcl_const1s / (ztcg * zfacx1s) + zterm1 = (zqsice[jk - 1, jl - 1] - + zqe) * ztp1[jk - 1, + jl - 1]**2 * zvpice * zcorrfac2 * ztcg * yrecldp_rcl_const2s * zfacx1s / ( + zrho[jl - 1] * zaplusb * zqsice[jk - 1, jl - 1]) + zterm2 = 0.65 * yrecldp_rcl_const6s * zpr02**yrecldp_rcl_const4s + yrecldp_rcl_const3s * zcorrfac**0.5 * zrho[ + jl - 1]**0.5 * zpr02**yrecldp_rcl_const5s / zcorrfac2**0.5 + zdpevap = max(zcovpclr[jl - 1] * zterm1 * zterm2 * ptsphy, 0.0) + zevap = min(zdpevap, zevaplimice[jl - 1]) + zevap = min(zevap, zqx[ncldqs - 1, jk - 1, jl - 1]) + zsolqa[ncldqs - 1, ncldqv - 1, jl - 1] = zsolqa[ncldqs - 1, ncldqv - 1, jl - 1] + zevap + zsolqa[ncldqv - 1, ncldqs - 1, jl - 1] = zsolqa[ncldqv - 1, ncldqs - 1, jl - 1] - zevap + zcovptot[jl - 1] = max( + yrecldp_rcovpmin, zcovptot[jl - 1] - + max(0.0, (zcovptot[jl - 1] - za[jk - 1, jl - 1]) * zevap / zqx[ncldqs - 1, jk - 1, jl - 1])) + zqxfg[ncldqs - 1, jl - 1] = zqxfg[ncldqs - 1, jl - 1] - zevap + for jm in range(1, nclv + 1): + if llfall[jm - 1]: + for jl in range(kidia, kfdia + 1): + if zqxfg[jm - 1, jl - 1] < yrecldp_rlmin: + zsolqa[jm - 1, ncldqv - 1, jl - 1] = zsolqa[jm - 1, ncldqv - 1, jl - 1] + zqxfg[jm - 1, jl - 1] + zsolqa[ncldqv - 1, jm - 1, jl - 1] = zsolqa[ncldqv - 1, jm - 1, jl - 1] - zqxfg[jm - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + zanew = (za[jk - 1, jl - 1] + zsolac[jl - 1]) / (1.0 + zsolab[jl - 1]) + zanew = min(zanew, 1.0) + if zanew < yrecldp_ramin: + zanew = 0.0 + zda[jl - 1] = zanew - zaorig[jk - 1, jl - 1] + zanewm1[jl - 1] = zanew + for jm in range(1, nclv + 1): + for jn in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + llindex3[jm - 1, jn - 1, jl - 1] = False + for jl in range(kidia, kfdia + 1): + zsinksum[jm - 1, jl - 1] = 0.0 + for jm in range(1, nclv + 1): + for jn in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zsinksum[jm - 1, jl - 1] = zsinksum[jm - 1, jl - 1] - zsolqa[jn - 1, jm - 1, jl - 1] + for jm in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zmax = max(zqx[jm - 1, jk - 1, jl - 1], zepsec) + zrat = max(zsinksum[jm - 1, jl - 1], zmax) + zratio[jm - 1, jl - 1] = zmax / zrat + for jm in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zsinksum[jm - 1, jl - 1] = 0.0 + for jm in range(1, nclv + 1): + for jl in range(1, klon + 1): + psum_solqa[jl - 1] = 0.0 + for jn in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + psum_solqa[jl - 1] = psum_solqa[jl - 1] + zsolqa[jn - 1, jm - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + zsinksum[jm - 1, jl - 1] = zsinksum[jm - 1, jl - 1] - psum_solqa[jl - 1] + for jl in range(kidia, kfdia + 1): + zmm = max(zqx[jm - 1, jk - 1, jl - 1], zepsec) + zrr = max(zsinksum[jm - 1, jl - 1], zmm) + zratio[jm - 1, jl - 1] = zmm / zrr + for jl in range(kidia, kfdia + 1): + zzratio = zratio[jm - 1, jl - 1] + for jn in range(1, nclv + 1): + if zsolqa[jn - 1, jm - 1, jl - 1] < 0.0: + zsolqa[jn - 1, jm - 1, jl - 1] = zsolqa[jn - 1, jm - 1, jl - 1] * zzratio + zsolqa[jm - 1, jn - 1, jl - 1] = zsolqa[jm - 1, jn - 1, jl - 1] * zzratio + for jm in range(1, nclv + 1): + for jn in range(1, nclv + 1): + if jn == jm: + for jl in range(kidia, kfdia + 1): + zqlhs[jm - 1, jn - 1, jl - 1] = 1.0 + zfallsink[jm - 1, jl - 1] + for jo in range(1, nclv + 1): + zqlhs[jm - 1, jn - 1, + jl - 1] = zqlhs[jm - 1, jn - 1, jl - 1] + zsolqb[jn - 1, jo - 1, jl - 1] + else: + for jl in range(kidia, kfdia + 1): + zqlhs[jm - 1, jn - 1, jl - 1] = -zsolqb[jm - 1, jn - 1, jl - 1] + for jm in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zexplicit = 0.0 + for jn in range(1, nclv + 1): + zexplicit = zexplicit + zsolqa[jn - 1, jm - 1, jl - 1] + zqxn[jm - 1, jl - 1] = zqx[jm - 1, jk - 1, jl - 1] + zexplicit + for jn in range(1, nclv - 1 + 1): + for jm in range(jn + 1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zqlhs[jn - 1, jm - 1, jl - 1] = zqlhs[jn - 1, jm - 1, jl - 1] / zqlhs[jn - 1, jn - 1, jl - 1] + for jn in range(1, nclv - 1 + 1): + for jm in range(jn + 1, nclv + 1): + for ik in range(jn + 1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zqlhs[ik - 1, jm - 1, + jl - 1] = zqlhs[ik - 1, jm - 1, + jl - 1] - zqlhs[jn - 1, jm - 1, jl - 1] * zqlhs[ik - 1, jn - 1, jl - 1] + for jn in range(2, nclv + 1): + for jm in range(1, jn - 1 + 1): + for jl in range(kidia, kfdia + 1): + zqxn[jn - 1, jl - 1] = zqxn[jn - 1, jl - 1] - zqlhs[jm - 1, jn - 1, jl - 1] * zqxn[jm - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + zqxn[nclv - 1, jl - 1] = zqxn[nclv - 1, jl - 1] / zqlhs[nclv - 1, nclv - 1, jl - 1] + for jn in range(nclv - 1, 1 + -1, -1): + for jm in range(jn + 1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zqxn[jn - 1, jl - 1] = zqxn[jn - 1, jl - 1] - zqlhs[jm - 1, jn - 1, jl - 1] * zqxn[jm - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + zqxn[jn - 1, jl - 1] = zqxn[jn - 1, jl - 1] / zqlhs[jn - 1, jn - 1, jl - 1] + for jn in range(1, nclv - 1 + 1): + for jl in range(kidia, kfdia + 1): + if zqxn[jn - 1, jl - 1] < zepsec: + zqxn[ncldqv - 1, jl - 1] = zqxn[ncldqv - 1, jl - 1] + zqxn[jn - 1, jl - 1] + zqxn[jn - 1, jl - 1] = 0.0 + for jm in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zqxnm1[jm - 1, jl - 1] = zqxn[jm - 1, jl - 1] + zqxn2d[jm - 1, jk - 1, jl - 1] = zqxn[jm - 1, jl - 1] + for jm in range(1, nclv + 1): + for jl in range(kidia, kfdia + 1): + zpfplsx[jm - 1, jk + 1 - 1, jl - 1] = zfallsink[jm - 1, jl - 1] * zqxn[jm - 1, jl - 1] * zrdtgdp[jl - 1] + for jl in range(kidia, kfdia + 1): + zqpretot[jl - 1] = zpfplsx[ncldqs - 1, jk + 1 - 1, jl - 1] + zpfplsx[ncldqr - 1, jk + 1 - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + if zqpretot[jl - 1] < zepsec: + zcovptot[jl - 1] = 0.0 + for jm in range(1, nclv - 1 + 1): + for jl in range(kidia, kfdia + 1): + zfluxq[jm - 1, jl - 1] = zpsupsatsrce[jm - 1, jl - 1] + zconvsrce[jm - 1, jl - 1] + zfallsrce[ + jm - 1, jl - 1] - (zfallsink[jm - 1, jl - 1] + zconvsink[jm - 1, jl - 1]) * zqxn[jm - 1, jl - 1] + if iphase[jm - 1] == 1: + for jl in range(kidia, kfdia + 1): + tendency_loc_t[jk - 1, jl - 1] = tendency_loc_t[jk - 1, jl - 1] + ydthf_ralvdcp * ( + zqxn[jm - 1, jl - 1] - zqx[jm - 1, jk - 1, jl - 1] - zfluxq[jm - 1, jl - 1]) * zqtmst + if iphase[jm - 1] == 2: + for jl in range(kidia, kfdia + 1): + tendency_loc_t[jk - 1, jl - 1] = tendency_loc_t[jk - 1, jl - 1] + ydthf_ralsdcp * ( + zqxn[jm - 1, jl - 1] - zqx[jm - 1, jk - 1, jl - 1] - zfluxq[jm - 1, jl - 1]) * zqtmst + for jl in range(kidia, kfdia + 1): + tendency_loc_cld[jm - 1, jk - 1, jl - + 1] = tendency_loc_cld[jm - 1, jk - 1, jl - 1] + (zqxn[jm - 1, jl - 1] - + zqx0[jm - 1, jk - 1, jl - 1]) * zqtmst + for jl in range(kidia, kfdia + 1): + tendency_loc_q[jk - 1, jl - 1] = tendency_loc_q[jk - 1, jl - 1] + (zqxn[ncldqv - 1, jl - 1] - + zqx[ncldqv - 1, jk - 1, jl - 1]) * zqtmst + tendency_loc_a[jk - 1, jl - 1] = tendency_loc_a[jk - 1, jl - 1] + zda[jl - 1] * zqtmst + for jl in range(kidia, kfdia + 1): + pcovptot[jk - 1, jl - 1] = zcovptot[jl - 1] + for jk in range(1, klev + 1 + 1): + for jl in range(kidia, kfdia + 1): + pfplsl[jk - 1, jl - 1] = zpfplsx[ncldqr - 1, jk - 1, jl - 1] + zpfplsx[ncldql - 1, jk - 1, jl - 1] + pfplsn[jk - 1, jl - 1] = zpfplsx[ncldqs - 1, jk - 1, jl - 1] + zpfplsx[ncldqi - 1, jk - 1, jl - 1] + for jl in range(kidia, kfdia + 1): + pfsqlf[1 - 1, jl - 1] = 0.0 + pfsqif[1 - 1, jl - 1] = 0.0 + pfsqrf[1 - 1, jl - 1] = 0.0 + pfsqsf[1 - 1, jl - 1] = 0.0 + pfcqlng[1 - 1, jl - 1] = 0.0 + pfcqnng[1 - 1, jl - 1] = 0.0 + pfcqrng[1 - 1, jl - 1] = 0.0 + pfcqsng[1 - 1, jl - 1] = 0.0 + pfsqltur[1 - 1, jl - 1] = 0.0 + pfsqitur[1 - 1, jl - 1] = 0.0 + for jk in range(1, klev + 1): + for jl in range(kidia, kfdia + 1): + zgdph_r = -zrg_r * (paph[jk + 1 - 1, jl - 1] - paph[jk - 1, jl - 1]) * zqtmst + pfsqlf[jk + 1 - 1, jl - 1] = pfsqlf[jk - 1, jl - 1] + pfsqif[jk + 1 - 1, jl - 1] = pfsqif[jk - 1, jl - 1] + pfsqrf[jk + 1 - 1, jl - 1] = pfsqlf[jk - 1, jl - 1] + pfsqsf[jk + 1 - 1, jl - 1] = pfsqif[jk - 1, jl - 1] + pfcqlng[jk + 1 - 1, jl - 1] = pfcqlng[jk - 1, jl - 1] + pfcqnng[jk + 1 - 1, jl - 1] = pfcqnng[jk - 1, jl - 1] + pfcqrng[jk + 1 - 1, jl - 1] = pfcqlng[jk - 1, jl - 1] + pfcqsng[jk + 1 - 1, jl - 1] = pfcqnng[jk - 1, jl - 1] + pfsqltur[jk + 1 - 1, jl - 1] = pfsqltur[jk - 1, jl - 1] + pfsqitur[jk + 1 - 1, jl - 1] = pfsqitur[jk - 1, jl - 1] + zalfaw = zfoealfa[jk - 1, jl - 1] + pfsqlf[jk + 1 - 1, + jl - 1] = pfsqlf[jk + 1 - 1, + jl - 1] + (zqxn2d[ncldql - 1, jk - 1, jl - 1] - zqx0[ncldql - 1, jk - 1, jl - 1] + + pvfl[jk - 1, jl - 1] * ptsphy - zalfaw * plude[jk - 1, jl - 1]) * zgdph_r + pfcqlng[jk + 1 - 1, jl - 1] = pfcqlng[jk + 1 - 1, jl - 1] + zlneg[ncldql - 1, jk - 1, jl - 1] * zgdph_r + pfsqltur[jk + 1 - 1, jl - 1] = pfsqltur[jk + 1 - 1, jl - 1] + pvfl[jk - 1, jl - 1] * ptsphy * zgdph_r + pfsqrf[jk + 1 - 1, jl - + 1] = pfsqrf[jk + 1 - 1, jl - + 1] + (zqxn2d[ncldqr - 1, jk - 1, jl - 1] - zqx0[ncldqr - 1, jk - 1, jl - 1]) * zgdph_r + pfcqrng[jk + 1 - 1, jl - 1] = pfcqrng[jk + 1 - 1, jl - 1] + zlneg[ncldqr - 1, jk - 1, jl - 1] * zgdph_r + pfsqif[jk + 1 - 1, jl - + 1] = pfsqif[jk + 1 - 1, jl - 1] + (zqxn2d[ncldqi - 1, jk - 1, jl - 1] - + zqx0[ncldqi - 1, jk - 1, jl - 1] + pvfi[jk - 1, jl - 1] * ptsphy - + (1.0 - zalfaw) * plude[jk - 1, jl - 1]) * zgdph_r + pfcqnng[jk + 1 - 1, jl - 1] = pfcqnng[jk + 1 - 1, jl - 1] + zlneg[ncldqi - 1, jk - 1, jl - 1] * zgdph_r + pfsqitur[jk + 1 - 1, jl - 1] = pfsqitur[jk + 1 - 1, jl - 1] + pvfi[jk - 1, jl - 1] * ptsphy * zgdph_r + pfsqsf[jk + 1 - 1, jl - + 1] = pfsqsf[jk + 1 - 1, jl - + 1] + (zqxn2d[ncldqs - 1, jk - 1, jl - 1] - zqx0[ncldqs - 1, jk - 1, jl - 1]) * zgdph_r + pfcqsng[jk + 1 - 1, jl - 1] = pfcqsng[jk + 1 - 1, jl - 1] + zlneg[ncldqs - 1, jk - 1, jl - 1] * zgdph_r + for jk in range(1, klev + 1 + 1): + for jl in range(kidia, kfdia + 1): + pfhpsl[jk - 1, jl - 1] = -ydcst_rlvtt * pfplsl[jk - 1, jl - 1] + pfhpsn[jk - 1, jl - 1] = -ydcst_rlstt * pfplsn[jk - 1, jl - 1] diff --git a/tests/corpus/cloudsc/generate_data_for_cloudsc.py b/tests/corpus/cloudsc/generate_data_for_cloudsc.py new file mode 100644 index 0000000000..d26c72d135 --- /dev/null +++ b/tests/corpus/cloudsc/generate_data_for_cloudsc.py @@ -0,0 +1,419 @@ +# Copyright 2019-2026 ETH Zurich and the DaCe authors. All rights reserved. +"""Input-data generation and run-and-compare helpers for the inlined CloudSC +kernel in ``tests/corpus/cloudsc/cloudsc.py``, for end-to-end numerical tests. + +The inlined ``cloudsc_py`` program requires no callbacks, so it compiles and +runs standalone. These helpers build a runnable SDFG, generate a +physically-realistic input set, and run two SDFGs on identical inputs to check +that a transformation (simplify, ConstantPropagation, ...) is numerically +faithful to a non-transformed reference. + +Typical use:: + + ref = build_cloudsc_sdfg(simplify=False) + sut = build_cloudsc_sdfg(simplify=False) + sut.simplify() + assert run_and_compare(ref, sut) + +Data generation **follows the dwarf-p-cloudsc reference dataset** +(``config-files/input.h5`` of the upstream ECMWF dwarf): the YDCST/YDTHF/YRECLDP +physical constants in :data:`CLOUDSC_CONSTANTS` are the exact values from that +file, and every input array is filled with random values inside the ``[min, +max]`` range observed there (:data:`CLOUDSC_INPUT_RANGES`), with the dwarf's +``ncldtop = 15`` so the vertical microphysics loop runs over a realistic depth. +We do **not** load ``input.h5`` at runtime -- the harness is self-contained and +only mirrors its values -- so it stays runnable without the external dataset. + +Filling thresholds and latent heats with the real constants (rather than +uniform ``[0, 1)`` noise) keeps the kernel out of its degenerate branch regimes, +which is what makes a transform-vs-reference comparison meaningful: a random +constant set sits on every ``MIN``/``MAX`` and ``< rlmin`` boundary, so a +harmless floating-point reassociation flips branches and masquerades as a bug. +""" +import copy +from typing import Dict, List, Optional, Tuple, Union + +import numpy as np + +import dace +from dace import dtypes, symbolic +from dace.sdfg import nodes +from tests.corpus.cloudsc.cloudsc import cloudsc_py + +#: Shape symbols and named integer index scalars. The cloud species ``ncldq*`` +#: and ``ncldtop`` match the dwarf-p-cloudsc reference; the grid is kept small +#: (``klev = klon = 32``) for a fast compiled run -- the physical input ranges +#: are bounds, not vertical profiles, so they stay valid at any grid size, and +#: ``ncldtop = 15`` still leaves a meaningful vertical microphysics loop. +CLOUDSC_SYMBOLS: Dict[str, int] = { + 'klev': 32, + 'klon': 32, + 'nclv': 5, + 'ncldql': 1, # liquid cloud water + 'ncldqi': 2, # ice cloud water + 'ncldqr': 3, # rain water + 'ncldqs': 4, # snow + 'ncldqv': 5, # vapour + 'kidia': 1, + 'kfdia': 32, +} + +#: Exact YDCST/YDTHF/YRECLDP constants from the dwarf-p-cloudsc ``input.h5`` +#: reference. The ``yrecldp_nssopt``/``ncldtop``/``laeri*`` entries are integer +#: scalars (cast on use); the rest are doubles. Mirrored here so the harness +#: needs no external dataset. +CLOUDSC_CONSTANTS: Dict[str, float] = { + 'ptsphy': 3600.0, + 'ydcst_rcpd': 1004.7088578330674, + 'ydcst_rd': 287.0596736665907, + 'ydcst_retv': 0.6077667316114637, + 'ydcst_rg': 9.80665, + 'ydcst_rlmlt': 333700.0, + 'ydcst_rlstt': 2834500.0, + 'ydcst_rlvtt': 2500800.0, + 'ydcst_rtt': 273.16, + 'ydcst_rv': 461.5249933083879, + 'ydthf_r2es': 380.1608703442847, + 'ydthf_r3ies': 22.587, + 'ydthf_r3les': 17.502, + 'ydthf_r4ies': -0.7, + 'ydthf_r4les': 32.19, + 'ydthf_r5alscp': 17451123.253362577, + 'ydthf_r5alvcp': 10497584.68169531, + 'ydthf_r5ies': 6185.67582, + 'ydthf_r5les': 4217.45694, + 'ydthf_ralfdcp': 332.1360187066693, + 'ydthf_ralsdcp': 2821.2152982440934, + 'ydthf_ralvdcp': 2489.0792795374246, + 'ydthf_rkoop1': 2.583, + 'ydthf_rkoop2': 0.0048116, + 'ydthf_rtice': 250.16000000000003, + 'ydthf_rticecu': 250.16000000000003, + 'ydthf_rtwat': 273.16, + 'ydthf_rtwat_rtice_r': 0.043478260869565216, + 'ydthf_rtwat_rticecu_r': 0.043478260869565216, + 'yrecldp_laericeauto': 0, + 'yrecldp_laericesed': 0, + 'yrecldp_laerliqautolsp': 0, + 'yrecldp_laerliqcoll': 0, + 'yrecldp_ncldtop': 15, + 'yrecldp_nssopt': 1, + 'yrecldp_ramid': 0.8, + 'yrecldp_ramin': 1e-08, + 'yrecldp_rccn': 125.0, + 'yrecldp_rcl_apb1': 714000000000.0, + 'yrecldp_rcl_apb2': 116000000.0, + 'yrecldp_rcl_apb3': 241.6, + 'yrecldp_rcl_cdenom1': 557000000000.0, + 'yrecldp_rcl_cdenom2': 103000000.0, + 'yrecldp_rcl_cdenom3': 204.0, + 'yrecldp_rcl_const1i': 3.6231880115136998e-06, + 'yrecldp_rcl_const1r': 1.382300767579509, + 'yrecldp_rcl_const1s': 3.6231880115136998e-06, + 'yrecldp_rcl_const2i': 6283185.307179586, + 'yrecldp_rcl_const2r': 2143.2299120517614, + 'yrecldp_rcl_const2s': 6283185.307179586, + 'yrecldp_rcl_const3i': 596.9998475835998, + 'yrecldp_rcl_const3r': 0.6349999999999998, + 'yrecldp_rcl_const3s': 596.9998475835998, + 'yrecldp_rcl_const4i': 0.6666666666666666, + 'yrecldp_rcl_const4r': -0.20000000000000018, + 'yrecldp_rcl_const4s': 0.6666666666666666, + 'yrecldp_rcl_const5i': 0.9211666666666667, + 'yrecldp_rcl_const5r': 8685252.965082133, + 'yrecldp_rcl_const5s': 0.9211666666666667, + 'yrecldp_rcl_const6i': 1.0000000948961185, + 'yrecldp_rcl_const6r': -4.8, + 'yrecldp_rcl_const6s': 1.0000000948961185, + 'yrecldp_rcl_const7s': 90363515.76351073, + 'yrecldp_rcl_const8s': 1.1756666666666666, + 'yrecldp_rcl_fac1': 4146.902789847063, + 'yrecldp_rcl_fac2': 0.5555555555555556, + 'yrecldp_rcl_fzrab': -0.66, + 'yrecldp_rcl_ka273': 0.024, + 'yrecldp_rcl_kk_cloud_num_land': 300.0, + 'yrecldp_rcl_kk_cloud_num_sea': 50.0, + 'yrecldp_rcl_kkaac': 67.0, + 'yrecldp_rcl_kkaau': 1350.0, + 'yrecldp_rcl_kkbac': 1.15, + 'yrecldp_rcl_kkbaun': -1.79, + 'yrecldp_rcl_kkbauq': 2.47, + 'yrecldp_rclcrit_land': 0.00055, + 'yrecldp_rclcrit_sea': 0.00025, + 'yrecldp_rcldiff': 3e-06, + 'yrecldp_rcldiff_convi': 7.0, + 'yrecldp_rcldtopcf': 0.01, + 'yrecldp_rcovpmin': 0.1, + 'yrecldp_rdensref': 1.0, + 'yrecldp_rdepliqrefdepth': 500.0, + 'yrecldp_rdepliqrefrate': 0.1, + 'yrecldp_riceinit': 1e-12, + 'yrecldp_rkconv': 0.00016666666666666666, + 'yrecldp_rkooptau': 10800.0, + 'yrecldp_rlcritsnow': 3e-05, + 'yrecldp_rlmin': 1e-08, + 'yrecldp_rnice': 0.027, + 'yrecldp_rpecons': 5.54725619859993e-05, + 'yrecldp_rprc1': 100.0, + 'yrecldp_rprecrhmax': 0.7, + 'yrecldp_rsnowlin1': 0.001, + 'yrecldp_rsnowlin2': 0.03, + 'yrecldp_rtaumel': 7200.0, + 'yrecldp_rthomo': 235.16000000000003, + 'yrecldp_rvice': 0.13, + 'yrecldp_rvrain': 4.0, + 'yrecldp_rvrfactor': 0.00509, + 'yrecldp_rvsnow': 1.0, +} + +#: ``[min, max]`` range of each floating input array in the dwarf reference. +#: Arrays absent here are kernel outputs (or have no reference value) and are +#: zero-initialized. Several reference inputs are uniformly zero (``pmfd``, +#: ``pnice``, ``plsm``, ...); their ``(0.0, 0.0)`` range reproduces that. +CLOUDSC_INPUT_RANGES: Dict[str, Tuple[float, float]] = { + 'pa': (0.0, 1.0), + 'pap': (0.999923895048429, 101254.97084602337), + 'paph': (0.0, 101375.10057152908), + 'pccn': (0.0, 0.0), + 'pclv': (0.0, 4.0253768484705866e-05), + 'pdyna': (-0.0001187964540362494, 0.00013518935141888), + 'pdyni': (-2.140507622298193e-09, 1.5815637844680652e-09), + 'pdynl': (-4.134462057260392e-09, 2.9599904813560976e-09), + 'phrlw': (-0.00016678085208516417, 1.6941956516047845e-05), + 'phrsw': (-4.112083689734362e-20, 0.0), + 'picrit_aer': (0.0, 0.0), + 'plcrit_aer': (0.0, 0.0), + 'plsm': (0.0, 0.0), + 'plu': (0.0, 0.00028126794898995864), + 'plude': (0.0, 1.2765809819686113e-06), + 'pmfd': (0.0, 0.0), + 'pmfu': (0.0, 0.06243987753452692), + 'pnice': (0.0, 0.0), + 'pq': (1.030691002421725e-06, 0.0024358218045027608), + 'pre_ice': (0.0, 0.0), + 'psnde': (0.0, 0.0), + 'psupsat': (0.0, 4.762341884143876e-05), + 'pt': (196.49936539855418, 267.57212905220615), + 'pvervel': (-0.20627060482285262, 0.1776496848080718), + 'pvfa': (-0.0002305417072000441, 0.0002777777777777778), + 'pvfi': (-3.6683428537982023e-09, 1.2988008242524594e-08), + 'pvfl': (-2.193240268281019e-09, 6.007184256012529e-09), + 'tendency_tmp_a': (-0.00024253423245031549, 0.0002777777777777778), + 'tendency_tmp_cld': (-4.164243445880186e-09, 1.2724998598908134e-08), + 'tendency_tmp_q': (-6.92205735885683e-08, 4.265980175425775e-08), + 'tendency_tmp_t': (-0.000472018266728236, 0.0003986018193989591), +} + +#: ``[min, max]`` integer range of each integer input array in the reference. +CLOUDSC_INT_RANGES: Dict[str, Tuple[int, int]] = { + 'ktype': (0, 3), + 'ldcum': (0, 1), +} + +#: Compiler flags for a deterministic IEEE build: ``-O0`` with no fast-math and +#: no FP contraction (FMA), so the compiler performs no floating-point +#: reassociation. Under these flags the simplified and reference cloudsc SDFGs +#: agree **bit-for-bit**, which is why the correctness check (see +#: :func:`run_and_compare`) uses this build with a ``1e-15`` tolerance. +IEEE_CPU_ARGS: str = '-std=c++14 -fPIC -O0 -fopenmp -fno-fast-math -ffp-contract=off' + +#: Compiler flags for an ``-O3`` build that still preserves IEEE semantics: +#: same ``-fno-fast-math -ffp-contract=off`` flags as the IEEE build, so the +#: compiler is free to schedule / unroll / vectorise but cannot reassociate +#: or fuse FP ops. This is the regime we actually want to validate cloudsc +#: under: fast-math (``-ffast-math``) lets the compiler rewrite transcendentals +#: and reorder reductions, which produces large drifts the test cannot bound +#: cleanly -- so the suite never enables it. +O3_CPU_ARGS: str = '-std=c++14 -fPIC -O3 -fopenmp -fno-fast-math -ffp-contract=off' + +PARSED_CLOUDSC: Dict[bool, dace.SDFG] = {} + + +def build_cloudsc_sdfg(simplify: bool = False) -> dace.SDFG: + """A private copy of the CloudSC SDFG: parsed once per process, deepcopied per caller. + + The deepcopy is not optional: the parse costs minutes so it is memoized, and every consumer + transforms the SDFG it gets. + """ + if simplify not in PARSED_CLOUDSC: + sdfg = cloudsc_py.to_sdfg(simplify=simplify) + sdfg.validate() + PARSED_CLOUDSC[simplify] = sdfg + return copy.deepcopy(PARSED_CLOUDSC[simplify]) + + +def generate_cloudsc_inputs(sdfg: dace.SDFG, seed: int = 0) -> Dict[str, Union[np.ndarray, int, float]]: + """Generate a physically-realistic CloudSC input set for ``sdfg``. + + Every non-transient argument is filled following the dwarf reference (see + the module docstring): named constants from :data:`CLOUDSC_CONSTANTS`, + floating arrays uniform within their :data:`CLOUDSC_INPUT_RANGES` window + (kernel outputs / unknown arrays zeroed), integer arrays uniform within + :data:`CLOUDSC_INT_RANGES`, and the named index scalars / shape symbols from + :data:`CLOUDSC_SYMBOLS`. Length-1 arrays are passed as scalars. + + :param sdfg: The CloudSC SDFG whose non-transient arrays are filled. + :param seed: Seed for the random number generator (reproducible runs). + :returns: A kwargs dict of arrays, scalars, and symbol values. + """ + rng = np.random.default_rng(seed) + arrays: Dict[str, np.ndarray] = {} + for name, desc in sdfg.arrays.items(): + if desc.transient: + continue + dims: List[int] = [int(symbolic.evaluate(d, CLOUDSC_SYMBOLS)) for d in desc.shape] + is_int = 'int' in str(desc.dtype) + + if name in CLOUDSC_CONSTANTS: + value = CLOUDSC_CONSTANTS[name] + arrays[name] = np.full(dims, + int(value) if is_int else value, + dtype=np.int32 if is_int else np.float64, + order='F') + elif is_int: + if name in CLOUDSC_SYMBOLS: + data = np.zeros(dims, order='F').astype(np.int32) + data.flat[0] = CLOUDSC_SYMBOLS[name] + arrays[name] = data + else: + lo, hi = CLOUDSC_INT_RANGES.get(name, (1, 1)) + arrays[name] = rng.integers(lo, hi + 1, size=dims).astype(np.int32, order='F') + else: + value_range: Optional[Tuple[float, float]] = CLOUDSC_INPUT_RANGES.get(name) + if value_range is None: + arrays[name] = np.zeros(dims, order='F') # kernel output / no reference range + else: + lo, hi = value_range + arrays[name] = (lo + (hi - lo) * rng.random(dims)).astype(np.float64, order='F') + + inputs: Dict[str, Union[np.ndarray, int, float]] = { + name: (data.flat[0] if data.size == 1 else data) + for name, data in arrays.items() + } + inputs.update(CLOUDSC_SYMBOLS) + return inputs + + +def make_sequential(sdfg: dace.SDFG) -> Tuple[int, int]: + """Force every map and library node to a sequential schedule, in place. + + A numerical-equivalence comparison must be deterministic: CloudSC's + parallel (OpenMP) maps reorder floating-point reductions / accumulations + run-to-run, so two separately-compiled SDFGs differ by ~1e-5 even when they + are the *same* computation. Running sequentially removes that noise so a + real numerical difference between a transform and its reference stands out. + + :param sdfg: The SDFG to make sequential (mutated in place). + :returns: ``(maps, library_nodes)`` re-scheduled. + """ + n_maps = n_lib = 0 + for node, _ in sdfg.all_nodes_recursive(): + if isinstance(node, nodes.MapEntry): + node.map.schedule = dtypes.ScheduleType.Sequential + n_maps += 1 + elif isinstance(node, nodes.LibraryNode): + node.schedule = dtypes.ScheduleType.Sequential + n_lib += 1 + return n_maps, n_lib + + +def compare_outputs(ref_inputs: Dict[str, Union[np.ndarray, int, float]], + cand_inputs: Dict[str, Union[np.ndarray, int, float]], + rtol: float = 1e-15, + atol: float = 1e-15, + verbose: bool = False) -> Dict[str, Tuple[float, float, bool]]: + """Per-array (per-subset) tolerance check between two driven input/output sets. + + Returns ``{array_name: (max_abs, max_rel, ok)}`` for every shared + ``numpy`` array, where ``max_abs`` is the worst absolute difference, + ``max_rel`` the worst relative difference over the nonzero reference + elements, and ``ok`` the per-array ``numpy.allclose`` verdict at ``rtol`` / + ``atol``. With ``verbose`` it prints one line per array. + + With these physical inputs and an IEEE build (:data:`IEEE_CPU_ARGS`), simplify + (and other value-preserving transforms) reproduce the no-transform reference + *bit-for-bit*, so the machine-precision ``1e-15`` default is appropriate for + a correctness check. A release build (``-O3 -ffast-math``, run in parallel) + approximates the transcendental intrinsics and reorders the flux prefix sums + and parallel reductions, so the caller should pass a looser tolerance + (~``1e-10``) there. + + :param ref_inputs: The reference run's mutated input/output dict. + :param cand_inputs: The candidate run's mutated input/output dict. + :param rtol: Per-array relative tolerance. + :param atol: Per-array absolute tolerance. + :param verbose: Print a ``name max_abs max_rel ok`` line per array. + :returns: Mapping from array name to ``(max_abs, max_rel, ok)``. + """ + report: Dict[str, Tuple[float, float, bool]] = {} + for name, ref_val in ref_inputs.items(): + if not isinstance(ref_val, np.ndarray) or ref_val.size == 0: + continue + cand_val = cand_inputs[name] + absd = np.abs(ref_val - cand_val) + max_abs = float(np.nanmax(absd)) + denom = np.abs(ref_val) + nz = denom > 0 + max_rel = float(np.nanmax(absd[nz] / denom[nz])) if np.any(nz) else 0.0 + ok = bool(np.allclose(ref_val, cand_val, rtol=rtol, atol=atol, equal_nan=True)) + report[name] = (max_abs, max_rel, ok) + if verbose: + print(f"{name:24s} max_abs={max_abs:.3e} max_rel={max_rel:.3e} {'ok' if ok else 'FAIL'}") + return report + + +def run_and_compare(reference: dace.SDFG, + candidate: dace.SDFG, + seed: int = 0, + rtol: float = 1e-15, + atol: float = 1e-15, + sequential: bool = True, + ieee_build: bool = True, + verbose: bool = False) -> bool: + """Run two CloudSC SDFGs on identical inputs and compare every shared + non-transient output array, each at its own ``rtol`` / ``atol``. + + Both SDFGs are driven with the same generated inputs (a private deep copy + each, since the call mutates the buffers in place). Only array outputs are + compared; scalar and symbol inputs are not. + + For a correctness check the default ``ieee_build`` compiles both SDFGs with + :data:`IEEE_CPU_ARGS` (``-O0``, no fast-math, no FP contraction), under which + a value-preserving transform reproduces the reference *bit-for-bit*, so the + default ``1e-15`` tolerance holds. To instead measure a release build, pass + ``ieee_build=False`` (leaving the configured ``-O3 -ffast-math`` flags) with + ``sequential=False`` and a looser tolerance (~``1e-10``): the transcendental + intrinsics are then approximated and the flux prefix sums and parallel + reductions reordered. + + :param reference: The reference SDFG (e.g. un-transformed). + :param candidate: The SDFG under test. + :param seed: Seed for input generation. + :param rtol: Per-array relative tolerance (``numpy.allclose``). + :param atol: Per-array absolute tolerance (``numpy.allclose``). + :param sequential: Force both SDFGs sequential before running + (:func:`make_sequential`) so the comparison is deterministic; leave on + unless you specifically want to measure parallel behaviour. + :param ieee_build: Compile with the deterministic IEEE flags + (:data:`IEEE_CPU_ARGS`) for the duration of the runs; the prior + ``compiler.cpu.args`` setting is restored afterwards. + :param verbose: Print the per-array ``max_abs`` / ``max_rel`` and verdict. + :returns: ``True`` iff every shared output array matches within tolerance. + """ + if sequential: + make_sequential(reference) + make_sequential(candidate) + ref_inputs = generate_cloudsc_inputs(reference, seed) + cand_inputs = copy.deepcopy(ref_inputs) + + saved_args = dace.Config.get('compiler', 'cpu', 'args') + try: + if ieee_build: + dace.Config.set('compiler', 'cpu', 'args', value=IEEE_CPU_ARGS) + reference(**ref_inputs) + candidate(**cand_inputs) + finally: + dace.Config.set('compiler', 'cpu', 'args', value=saved_args) + + report = compare_outputs(ref_inputs, cand_inputs, rtol=rtol, atol=atol, verbose=verbose) + return all(ok for _, _, ok in report.values()) diff --git a/tests/corpus/cloudsc_regression_test.py b/tests/corpus/cloudsc_regression_test.py new file mode 100644 index 0000000000..4bd24adf2a --- /dev/null +++ b/tests/corpus/cloudsc_regression_test.py @@ -0,0 +1,51 @@ +# Copyright 2019-2026 ETH Zurich and the DaCe authors. All rights reserved. +"""Wall-clock regression guard on simplifying the CloudSC corpus kernel. + +CloudSC is the scaling case: thousands of blocks nested many levels deep, where a pass whose cost +is superlinear in nesting depth shows up as minutes rather than seconds. ConstantPropagation is the +dominant term in ``simplify`` on this input. +""" +import statistics +import time + +import pytest + +import dace +from tests.corpus.cloudsc.generate_data_for_cloudsc import build_cloudsc_sdfg + +#: Wall-clock budget for one ``simplify()`` of CloudSC. +SIMPLIFY_BUDGET_SECONDS: float = 140.0 +SIMPLIFY_REPS: int = 5 + + +@pytest.mark.long +def test_build_cloudsc_sdfg_hands_out_private_copies(): + """The corpus memoizes the parse; every caller must still get a copy it can transform.""" + first = build_cloudsc_sdfg(simplify=False) + second = build_cloudsc_sdfg(simplify=False) + assert first is not second + first.add_symbol('canary', dace.int32) + assert 'canary' not in second.symbols + + +@pytest.mark.long +def test_simplify_stays_within_its_time_budget(): + """One run per copy: ``simplify`` is idempotent, so a second run on the same SDFG measures + nothing. The median absorbs a slow rep from a loaded machine without needing a wide margin. + """ + durations = [] + for _ in range(SIMPLIFY_REPS): + sdfg = build_cloudsc_sdfg(simplify=False) # deliberately outside the timer + start = time.perf_counter() + sdfg.simplify() + durations.append(time.perf_counter() - start) + + median = statistics.median(durations) + reps = ', '.join('%.1f' % d for d in durations) + assert median < SIMPLIFY_BUDGET_SECONDS, (f'median simplify took {median:.1f}s, budget is ' + f'{SIMPLIFY_BUDGET_SECONDS:.0f}s; reps={reps}') + + +if __name__ == '__main__': + test_build_cloudsc_sdfg_hands_out_private_copies() + test_simplify_stays_within_its_time_budget() diff --git a/tests/passes/constant_propagation_on_cloudsc_test.py b/tests/passes/constant_propagation_on_cloudsc_test.py new file mode 100644 index 0000000000..426ccd2e1e --- /dev/null +++ b/tests/passes/constant_propagation_on_cloudsc_test.py @@ -0,0 +1,73 @@ +# Copyright 2019-2026 ETH Zurich and the DaCe authors. All rights reserved. +"""ConstantPropagation on the CloudSC corpus kernel. + +CloudSC is the scaling case for this pass: thousands of blocks and loop regions nested many levels +deep, which is where the collection schedule dominates the runtime. The last test closes the loop +by compiling and running, so a schedule change that quietly altered a value cannot pass. +""" +import pytest + +import dace +from dace.transformation.passes.constant_propagation import ConstantPropagation +from dace.transformation.passes.scalar_to_symbol import ScalarToSymbolPromotion +from dace.transformation.passes.simplification.control_flow_raising import ControlFlowRaising +from tests.corpus.cloudsc.generate_data_for_cloudsc import build_cloudsc_sdfg, run_and_compare + + +def promoted_cloudsc() -> dace.SDFG: + """CloudSC as SimplifyPass hands it to ConstantPropagation. + + Without the promotion no symbol carries a constant and the pass early-exits. + """ + sdfg = build_cloudsc_sdfg(simplify=False) + ScalarToSymbolPromotion().apply_pass(sdfg, {}) + ControlFlowRaising().apply_pass(sdfg, {}) + return sdfg + + +@pytest.mark.long +def test_constant_propagation_cloudsc(): + sdfg = promoted_cloudsc() + + propagated = ConstantPropagation().apply_pass(sdfg, {}) + assert propagated, 'expected constants to propagate in cloudsc' + sdfg.validate() + + # No constant assignment may survive on an interstate edge for a symbol that was propagated away. + propagated_symbols = {sym for cfg_id, sym in propagated if cfg_id == sdfg.cfg_id} + for edge in sdfg.all_interstate_edges(): + assert not (propagated_symbols & edge.data.assignments.keys()) + + # An empty set here instead of None would spin SimplifyPass's FixedPointPipeline forever. + assert ConstantPropagation().apply_pass(sdfg, {}) is None + + +@pytest.mark.long +def test_constant_propagation_cloudsc_is_numerically_faithful(): + """End to end: promote, propagate, compile, run, compare against the un-propagated kernel.""" + reference = build_cloudsc_sdfg(simplify=False) + candidate = promoted_cloudsc() + # Under the 'name' cache config the build folder is just the SDFG name, so equal names collide. + candidate.name = f'{candidate.name}_propagated' + + assert ConstantPropagation().apply_pass(candidate, {}) + assert run_and_compare(reference, candidate) + + +@pytest.mark.long +def test_simplified_cloudsc_is_numerically_faithful(): + """The same check through the full ``simplify``, which runs ConstantPropagation in a pipeline.""" + reference = build_cloudsc_sdfg(simplify=False) + candidate = build_cloudsc_sdfg(simplify=False) + candidate.name = f'{candidate.name}_simplified' + + candidate.simplify(validate=True) + assert sum(1 for _ in candidate.all_control_flow_blocks()) < sum(1 for _ in reference.all_control_flow_blocks()) + + assert run_and_compare(reference, candidate) + + +if __name__ == '__main__': + test_constant_propagation_cloudsc() + test_constant_propagation_cloudsc_is_numerically_faithful() + test_simplified_cloudsc_is_numerically_faithful()