From 14af763dedb54b46657fd77b32f10403ae563ebf Mon Sep 17 00:00:00 2001 From: "ngo-nghi-truyen.huynh" Date: Wed, 31 Jul 2024 23:32:26 +0200 Subject: [PATCH] Fix PR from FC first review --- smash/_constant.py | 4 +- smash/core/model/_build_model.py | 4 +- smash/core/model/_standardize.py | 10 +- smash/core/model/model.py | 18 +- smash/core/simulation/_doc.py | 27 +- smash/core/simulation/_standardize.py | 27 +- .../core/simulation/optimize/_standardize.py | 6 - smash/factory/net/_loss.py | 4 +- smash/factory/net/net.py | 4 +- .../derived_type/mwd_optimize_options.f90 | 3 +- smash/fcore/forward/forward_db.f90 | 1409 ++++++---------- smash/fcore/forward/forward_openmp_db.f90 | 1435 ++++++----------- smash/fcore/operator/md_gr_operator.f90 | 167 +- 13 files changed, 1144 insertions(+), 1974 deletions(-) diff --git a/smash/_constant.py b/smash/_constant.py index 68b15812..23b53053 100644 --- a/smash/_constant.py +++ b/smash/_constant.py @@ -589,7 +589,7 @@ def get_rr_internal_fluxes_from_structure(structure: str) -> list[str]: ### OPTIMIZABLE NN PARAMETERS ### ################################# -OPTIMIZABLE_NN_PARAMETERS = ["weight_1", "bias_1", "weight_2", "bias_2"] +NN_PARAMETERS_KEYS = ["weight_1", "bias_1", "weight_2", "bias_2"] ### SETUP ### ############# @@ -995,7 +995,7 @@ def get_rr_internal_fluxes_from_structure(structure: str) -> list[str]: "atmos_data": ["mean_prcp", "mean_pet", "mean_snow", "mean_temp"], "rr_parameters": ["keys", "values"], "rr_initial_states": ["keys", "values"], - "nn_parameters": OPTIMIZABLE_NN_PARAMETERS, + "nn_parameters": NN_PARAMETERS_KEYS, "serr_mu_parameters": ["keys", "values"], "serr_sigma_parameters": ["keys", "values"], "response": ["q"], diff --git a/smash/core/model/_build_model.py b/smash/core/model/_build_model.py index 57dbe1b7..d5469704 100644 --- a/smash/core/model/_build_model.py +++ b/smash/core/model/_build_model.py @@ -12,7 +12,7 @@ DEFAULT_RR_PARAMETERS, DEFAULT_SERR_MU_PARAMETERS, DEFAULT_SERR_SIGMA_PARAMETERS, - OPTIMIZABLE_NN_PARAMETERS, + NN_PARAMETERS_KEYS, SERR_MU_MAPPING_PARAMETERS, SERR_SIGMA_MAPPING_PARAMETERS, STRUCTURE_ADJUST_CI, @@ -181,7 +181,7 @@ def _build_parameters( # % Initalize weights and biases of ANN if hybrid model structure is used if sum(setup.neurons) > 0: - for key in OPTIMIZABLE_NN_PARAMETERS: + for key in NN_PARAMETERS_KEYS: # zero init setattr(parameters.nn_parameters, key, 0) diff --git a/smash/core/model/_standardize.py b/smash/core/model/_standardize.py index 52efa62c..b5e7ca87 100644 --- a/smash/core/model/_standardize.py +++ b/smash/core/model/_standardize.py @@ -416,7 +416,7 @@ def _standardize_model_setup_descriptor_name(descriptor_name: ListLike | None, * def _standardize_model_setup_hidden_neuron(hidden_neuron: Numeric, **kwargs) -> int: - if isinstance(hidden_neuron, (int, float, np.number)): + if isinstance(hidden_neuron, (int, float)): hidden_neuron = int(hidden_neuron) if hidden_neuron <= 0: raise ValueError("hidden_neuron model setup must be a positive number") @@ -763,7 +763,9 @@ def _standardize_set_nn_parameters_weight_value( f"from shape {arr.shape} into shape {weights[i].shape}" ) else: - raise TypeError("Each element of value argument must be float or a Numpy array") + raise TypeError( + "Each element of value argument must be of Numeric type (int, float) or np.ndarray" + ) else: raise TypeError("value argument must be a list of a same size with layers") @@ -797,7 +799,9 @@ def _standardize_set_nn_parameters_bias_value( f"from shape {arr.shape} into shape {biases[i].shape}" ) else: - raise TypeError("Each element of value argument must be float or a Numpy array") + raise TypeError( + "Each element of value argument must be of Numeric type (int, float) or np.ndarray" + ) else: raise TypeError("value argument must be a list of a same size with layers") diff --git a/smash/core/model/model.py b/smash/core/model/model.py index bdfcfc39..1e43b93a 100644 --- a/smash/core/model/model.py +++ b/smash/core/model/model.py @@ -1145,7 +1145,7 @@ def nn_parameters(self) -> NN_ParametersDT: >>> setup["hydrological_module"] = "gr4_mlp_alg" >>> model = smash.Model(setup, mesh) - By default, the weights and biases of the parameterization neural network is set to zero. + By default, the weight and bias of the parameterization neural network are set to zero. Access to their values with the getter method `get_nn_parameters_weight ` or `get_nn_parameters_bias ` @@ -2394,7 +2394,7 @@ def get_nn_parameters_weight(self) -> list[NDArray[np.float32]]: Returns ------- value : list[`numpy.ndarray`] - A list of arrays representing the weights for each layer of the parameterization neural network. + A list of arrays representing the weights of all layers. See Also -------- @@ -2416,7 +2416,7 @@ def get_nn_parameters_weight(self) -> list[NDArray[np.float32]]: >>> setup["hidden_neuron"] = 3 >>> model = smash.Model(setup, mesh) - By default, the weights of the parameterization neural network is set to zero. + By default, the weights of all layers are set to zero. Access to their values with the getter methods `get_nn_parameters_weight ` @@ -2440,7 +2440,7 @@ def get_nn_parameters_bias(self) -> list[NDArray[np.float32]]: Returns ------- value : list[`numpy.ndarray`] - A list of arrays representing the biases for each layer of the parameterization neural network. + A list of arrays representing the biases of all layers. See Also -------- @@ -2462,7 +2462,7 @@ def get_nn_parameters_bias(self) -> list[NDArray[np.float32]]: >>> setup["hidden_neuron"] = 6 >>> model = smash.Model(setup, mesh) - By default, the biases of the parameterization neural network is set to zero. + By default, the biases of all layers are set to zero. Access to their values with the getter methods `get_nn_parameters_bias ` @@ -2486,8 +2486,8 @@ def set_nn_parameters_weight( Parameters ---------- value : list[`float` or `numpy.ndarray`] or None, default None - The list of value(s) to set to the weights of the neural network. - If the value is a `numpy.ndarray`, its shape must be broadcastable into the weight shape. + The list of values to set to the weights of all layers. If an element of the list is + a `numpy.ndarray`, its shape must be broadcastable into the weight shape of that layer. If not used, a default or specified initialization method will be used. initializer : str, default 'glorot_uniform' @@ -2580,8 +2580,8 @@ def set_nn_parameters_bias( Parameters ---------- value : list[`float` or `numpy.ndarray`] or None, default None - The list of value(s) to set to the biases of the neural network. - If the value is a `numpy.ndarray`, its shape must be broadcastable into the bias shape. + The list of values to set to the biases of all layers. If an element of the list is + a `numpy.ndarray`, its shape must be broadcastable into the bias shape of that layer. If not used, a default or specified initialization method will be used. initializer : str, default 'zeros' diff --git a/smash/core/simulation/_doc.py b/smash/core/simulation/_doc.py index 5dac2830..aa34ce7a 100644 --- a/smash/core/simulation/_doc.py +++ b/smash/core/simulation/_doc.py @@ -1434,7 +1434,11 @@ def _gen_docstring_from_base_doc( ``'llr-slope-b'``, ``'ht-dd-a'``). - ``-``: Structural error parameter where ```` is the name of any structural error mu or sigma parameters and ````, the corresponding gauge (``'sg0-V3524010'``, ``'sg1-V3524010'``, - etc) + etc). + - ``--``: Weights and biases of the parameterization neural network where ```` + indicates the layer and type of parameter (e.g., ``'weight_1'`` for the first layer weights, + ``'bias_2'`` for the second layer biases), and ````, ```` represent the corresponding + position in the matrix or vector (``'weight_2-23-21'``, ``'bias_1-16'``, etc). - x_bkg : `numpy.ndarray` An array of shape *(n,)* containing the background values of the control vector. @@ -1605,11 +1609,12 @@ def _gen_docstring_from_base_doc( _smash_bayesian_optimize_doc_substitution = DocSubstitution( model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model `smash`.", default_optimize_options_func="default_bayesian_optimize_options", - parameters_nn_parameters="", - parameters_note_nn_parameters="", + parameters_nn_parameters="- `Model.nn_parameters`, if using a hybrid structure model " + "(depending on **hydrological_module**)", + parameters_note_nn_parameters=", `Model.nn_parameters` (if used)", parameters_serr_mu_parameters="- `Model.serr_mu_parameters`", parameters_serr_sigma_parameters="- `Model.serr_sigma_parameters`", - parameters_note_serr_parameters=", `Model.serr_mu_parameters` and `Model.serr_sigma_parameters`", + parameters_note_serr_parameters=", `Model.serr_mu_parameters`, `Model.serr_sigma_parameters`", bounds_get_serr_parameters_bounds=", `Model.get_serr_mu_parameters_bounds` and " "`Model.get_serr_sigma_parameters_bounds`", model_return="model : `Model`\n\t It returns an updated copy of the initial Model object.", @@ -1620,11 +1625,12 @@ def _gen_docstring_from_base_doc( _model_bayesian_optimize_doc_substitution = DocSubstitution( model_parameter="", default_optimize_options_func="default_bayesian_optimize_options", - parameters_nn_parameters="", - parameters_note_nn_parameters="", + parameters_nn_parameters="- `Model.nn_parameters`, if using a hybrid structure model " + "(depending on **hydrological_module**)", + parameters_note_nn_parameters=", `Model.nn_parameters` (if used)", parameters_serr_mu_parameters="- `Model.serr_mu_parameters`", parameters_serr_sigma_parameters="- `Model.serr_sigma_parameters`", - parameters_note_serr_parameters=", `Model.serr_mu_parameters` and `Model.serr_sigma_parameters`", + parameters_note_serr_parameters=", `Model.serr_mu_parameters`, `Model.serr_sigma_parameters`", bounds_get_serr_parameters_bounds=", `Model.get_serr_mu_parameters_bounds` and " "`Model.get_serr_sigma_parameters_bounds`", model_return="", @@ -1650,11 +1656,12 @@ def _gen_docstring_from_base_doc( _bayesian_optimize_control_info_doc_appender = DocAppender(_bayesian_optimize_control_info_doc, indents=0) _smash_bayesian_optimize_control_info_doc_substitution = DocSubstitution( default_optimize_options_func="default_bayesian_optimize_options", - parameters_nn_parameters="", - parameters_note_nn_parameters="", + parameters_nn_parameters="- `Model.nn_parameters`, if using a hybrid structure model " + "(depending on **hydrological_module**)", + parameters_note_nn_parameters=", `Model.nn_parameters` (if used)", parameters_serr_mu_parameters="- `Model.serr_mu_parameters`", parameters_serr_sigma_parameters="- `Model.serr_sigma_parameters`", - parameters_note_serr_parameters=", `Model.serr_mu_parameters` and `Model.serr_sigma_parameters`", + parameters_note_serr_parameters=", `Model.serr_mu_parameters`, `Model.serr_sigma_parameters`", bounds_get_serr_parameters_bounds=", `Model.get_serr_mu_parameters_bounds` and " "`Model.get_serr_sigma_parameters_bounds`", ) diff --git a/smash/core/simulation/_standardize.py b/smash/core/simulation/_standardize.py index ab6dff55..9ade8134 100644 --- a/smash/core/simulation/_standardize.py +++ b/smash/core/simulation/_standardize.py @@ -27,7 +27,7 @@ MAPPING, MAPPING_OPTIMIZER, METRICS, - OPTIMIZABLE_NN_PARAMETERS, + NN_PARAMETERS_KEYS, OPTIMIZABLE_RR_INITIAL_STATES, OPTIMIZABLE_RR_PARAMETERS, OPTIMIZABLE_SERR_MU_PARAMETERS, @@ -146,24 +146,21 @@ def _standardize_simulation_optimize_options_parameters( available_parameters = available_rr_parameters + available_rr_initial_states if is_hybrid_structure: - available_parameters.extend(OPTIMIZABLE_NN_PARAMETERS) + available_parameters.extend(NN_PARAMETERS_KEYS) if is_bayesian: available_parameters.extend(available_serr_mu_parameters + available_serr_sigma_parameters) if parameters is None: + default_parameters = available_rr_parameters + if is_bayesian: - parameters = np.array( - available_rr_parameters + available_serr_mu_parameters + available_serr_sigma_parameters, - ndmin=1, - ) - elif ( - is_hybrid_structure - ): # hybrid structure is currently not working with bayes optim. TODO for later versions - parameters = np.array(available_rr_parameters + OPTIMIZABLE_NN_PARAMETERS, ndmin=1) + default_parameters += available_serr_mu_parameters + available_serr_sigma_parameters - else: - parameters = np.array(available_rr_parameters, ndmin=1) + if is_hybrid_structure: + default_parameters += NN_PARAMETERS_KEYS + + parameters = np.array(default_parameters, ndmin=1) else: if isinstance(parameters, (str, list, tuple, np.ndarray)): @@ -192,7 +189,7 @@ def _standardize_simulation_optimize_options_parameters( def _standardize_simulation_optimize_options_bounds( model: Model, parameters: np.ndarray, bounds: dict | None, **kwargs ) -> dict: - bounded_parameters = [p for p in parameters if p not in OPTIMIZABLE_NN_PARAMETERS] + bounded_parameters = [p for p in parameters if p not in NN_PARAMETERS_KEYS] if bounds is None: bounds = {} @@ -1207,9 +1204,9 @@ def _standardize_simulation_optimize_options_finalize( optimize_options["rr_initial_states_descriptor"][j, i] = 1 # % nn parameters - optimize_options["nn_parameters"] = np.zeros(shape=len(OPTIMIZABLE_NN_PARAMETERS), dtype=np.int32) + optimize_options["nn_parameters"] = np.zeros(shape=len(NN_PARAMETERS_KEYS), dtype=np.int32) - for i, key in enumerate(OPTIMIZABLE_NN_PARAMETERS): + for i, key in enumerate(NN_PARAMETERS_KEYS): if key in optimize_options["parameters"]: optimize_options["nn_parameters"][i] = 1 diff --git a/smash/core/simulation/optimize/_standardize.py b/smash/core/simulation/optimize/_standardize.py index f48f54b8..8231bbd8 100644 --- a/smash/core/simulation/optimize/_standardize.py +++ b/smash/core/simulation/optimize/_standardize.py @@ -101,12 +101,6 @@ def _standardize_bayesian_optimize_args( ) -> AnyTuple: func_name = "bayesian_optimize" - if "mlp" in model.setup.hydrological_module: - raise ValueError( - f"Bayesian optimization method is currently not working with the hybrid model structure " - f"'{model.setup.hydrological_module}'" - ) - # % In case model.set_rr_parameters or model.set_rr_initial_states were not used _standardize_simulation_parameters_feasibility(model) diff --git a/smash/factory/net/_loss.py b/smash/factory/net/_loss.py index 02db4315..6d549e2a 100644 --- a/smash/factory/net/_loss.py +++ b/smash/factory/net/_loss.py @@ -4,7 +4,7 @@ import numpy as np -from smash._constant import OPTIMIZABLE_NN_PARAMETERS +from smash._constant import NN_PARAMETERS_KEYS from smash.fcore._mw_forward import forward_run_b as wrap_forward_run_b from smash.fcore._mwd_parameters_manipulation import ( control_to_parameters as wrap_control_to_parameters, @@ -102,7 +102,7 @@ def _hcost_prime( grad_reg = np.reshape(grad_reg, (len(grad_reg), -1)).T # % Get the gradient of parameterization NN if used - grad_par = [getattr(parameters_b.nn_parameters, key).copy() for key in OPTIMIZABLE_NN_PARAMETERS] + grad_par = [getattr(parameters_b.nn_parameters, key).copy() for key in NN_PARAMETERS_KEYS] return grad_reg, grad_par diff --git a/smash/factory/net/net.py b/smash/factory/net/net.py index c79451b4..85827c08 100644 --- a/smash/factory/net/net.py +++ b/smash/factory/net/net.py @@ -2,7 +2,7 @@ from typing import TYPE_CHECKING -from smash._constant import OPTIMIZABLE_NN_PARAMETERS, PY_OPTIMIZER, PY_OPTIMIZER_CLASS +from smash._constant import NN_PARAMETERS_KEYS, PY_OPTIMIZER, PY_OPTIMIZER_CLASS from smash.factory.net._layers import Activation, Conv2D, Dense, Dropout, Flatten, Scale from smash.factory.net._loss import _hcost, _hcost_prime, _inf_norm @@ -564,7 +564,7 @@ def _fit_d2p( # backpropagation and weights update if epo < epochs - 1: - for i, key in enumerate(OPTIMIZABLE_NN_PARAMETERS): + for i, key in enumerate(NN_PARAMETERS_KEYS): if key in parameters: # update trainable parameters of the parameterization NN if used setattr( instance.nn_parameters, diff --git a/smash/fcore/derived_type/mwd_optimize_options.f90 b/smash/fcore/derived_type/mwd_optimize_options.f90 index e8dc1f0c..6d3ae6b7 100644 --- a/smash/fcore/derived_type/mwd_optimize_options.f90 +++ b/smash/fcore/derived_type/mwd_optimize_options.f90 @@ -99,8 +99,7 @@ subroutine Optimize_OptionsDT_initialise(this, setup) allocate (this%rr_parameters_descriptor(setup%nd, setup%nrrp)) this%rr_parameters_descriptor = -99 - !% contain 4 matrices/vectors: weight_1, bias_1, weight_2, bias_2 - allocate (this%nn_parameters(4)) + allocate (this%nn_parameters((size(setup%neurons) - 1)*2)) this%nn_parameters = -99 allocate (this%rr_initial_states(setup%nrrs)) diff --git a/smash/fcore/forward/forward_db.f90 b/smash/fcore/forward/forward_db.f90 index ed151dd4..f7c128e2 100644 --- a/smash/fcore/forward/forward_db.f90 +++ b/smash/fcore/forward/forward_db.f90 @@ -12585,7 +12585,7 @@ SUBROUTINE DOT_PRODUCT_2D_1D_D(a, a_d, x, x_d, b, b_d) END SUBROUTINE DOT_PRODUCT_2D_1D_D ! Differentiation of dot_product_2d_1d in reverse (adjoint) mode (with options fixinterface noISIZE context): -! gradient of useful results: x a b +! gradient of useful results: a b ! with respect to varying inputs: x a SUBROUTINE DOT_PRODUCT_2D_1D_B(a, a_b, x, x_b, b, b_b) IMPLICIT NONE @@ -12606,6 +12606,7 @@ SUBROUTINE DOT_PRODUCT_2D_1D_B(a, a_b, x, x_b, b, b_b) CALL PUSHINTEGER4(i - 1) END DO ad_to0 = j - 1 + x_b = 0.0_4 DO j=ad_to0,1,-1 CALL POPINTEGER4(ad_to) DO i=ad_to,1,-1 @@ -12696,7 +12697,7 @@ END SUBROUTINE FORWARD_MLP_D ! Differentiation of forward_mlp in reverse (adjoint) mode (with options fixinterface noISIZE context): ! gradient of useful results: output_layer bias_1 bias_2 -! input_layer weight_1 weight_2 +! weight_1 weight_2 ! with respect to varying inputs: bias_1 bias_2 input_layer weight_1 ! weight_2 SUBROUTINE FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & @@ -12748,7 +12749,6 @@ SUBROUTINE FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & output_layer_b(i) = temp_b bias_2_b(i) = bias_2_b(i) + temp_b END DO - inter_layer_b = 0.0_4 CALL DOT_PRODUCT_2D_1D_B(weight_2, weight_2_b, inter_layer, & & inter_layer_b, output_layer, output_layer_b) CALL POPINTEGER4(ad_to) @@ -12803,7 +12803,6 @@ END MODULE MD_NEURAL_NETWORK_DIFF !% - gr_exchange !% - gr_threshold_exchange !% - gr_transfer -!% - gr_production_transfer_mlp_alg !% - gr_production_transfer_ode !% - gr_production_transfer_mlp_ode !% - gr4_time_step @@ -12955,12 +12954,12 @@ END SUBROUTINE GR_INTERCEPTION ! Differentiation of gr_production in forward (tangent) mode (with options fixinterface noISIZE context): ! variations of useful results: hp perc pr -! with respect to varying inputs: hp en cp pn - SUBROUTINE GR_PRODUCTION_D(pn, pn_d, en, en_d, cp, cp_d, beta, hp, & -& hp_d, pr, pr_d, perc, perc_d) +! with respect to varying inputs: fq_ps hp en fq_es cp pn + SUBROUTINE GR_PRODUCTION_D(fq_ps, fq_ps_d, fq_es, fq_es_d, pn, pn_d, & +& en, en_d, cp, cp_d, beta, hp, hp_d, pr, pr_d, perc, perc_d) IMPLICIT NONE - REAL(sp), INTENT(IN) :: pn, en, cp, beta - REAL(sp), INTENT(IN) :: pn_d, en_d, cp_d + REAL(sp), INTENT(IN) :: fq_ps, fq_es, pn, en, cp, beta + REAL(sp), INTENT(IN) :: fq_ps_d, fq_es_d, pn_d, en_d, cp_d REAL(sp), INTENT(INOUT) :: hp REAL(sp), INTENT(INOUT) :: hp_d REAL(sp), INTENT(OUT) :: pr, perc @@ -12987,6 +12986,8 @@ SUBROUTINE GR_PRODUCTION_D(pn, pn_d, en, en_d, cp, cp_d, beta, hp, & & inv_cp)**2)*(inv_cp*pn_d+pn*inv_cp_d)-temp2*(temp*hp_d+hp*(1.0-& & TANH(pn*inv_cp)**2)*(inv_cp*pn_d+pn*inv_cp_d)))/(hp*temp+1._sp) ps = temp2 + ps_d = ps*fq_ps_d + (fq_ps+1._sp)*ps_d + ps = (1._sp+fq_ps)*ps temp2 = TANH(en*inv_cp) temp1 = TANH(en*inv_cp) temp0 = hp*cp*(-hp+2._sp) @@ -12996,6 +12997,8 @@ SUBROUTINE GR_PRODUCTION_D(pn, pn_d, en, en_d, cp, cp_d, beta, hp, & & 1.0-TANH(en*inv_cp)**2)*(inv_cp*en_d+en*inv_cp_d)-temp2*hp_d))/((& & 1._sp-hp)*temp2+1._sp) es = temp + es_d = es*fq_es_d + (fq_es+1._sp)*es_d + es = (1._sp+fq_es)*es hp_imd_d = hp_d + inv_cp*(ps_d-es_d) + (ps-es)*inv_cp_d hp_imd = hp + (ps-es)*inv_cp IF (pn .GT. 0) THEN @@ -13015,13 +13018,14 @@ SUBROUTINE GR_PRODUCTION_D(pn, pn_d, en, en_d, cp, cp_d, beta, hp, & END SUBROUTINE GR_PRODUCTION_D ! Differentiation of gr_production in reverse (adjoint) mode (with options fixinterface noISIZE context): -! gradient of useful results: hp cp perc pr -! with respect to varying inputs: hp en cp pn - SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & -& hp_b, pr, pr_b, perc, perc_b) - IMPLICIT NONE - REAL(sp), INTENT(IN) :: pn, en, cp, beta - REAL(sp) :: pn_b, en_b, cp_b +! gradient of useful results: fq_ps hp en fq_es cp pn perc +! pr +! with respect to varying inputs: fq_ps hp en fq_es cp pn + SUBROUTINE GR_PRODUCTION_B(fq_ps, fq_ps_b, fq_es, fq_es_b, pn, pn_b, & +& en, en_b, cp, cp_b, beta, hp, hp_b, pr, pr_b, perc, perc_b) + IMPLICIT NONE + REAL(sp), INTENT(IN) :: fq_ps, fq_es, pn, en, cp, beta + REAL(sp) :: fq_ps_b, fq_es_b, pn_b, en_b, cp_b REAL(sp), INTENT(INOUT) :: hp REAL(sp), INTENT(INOUT) :: hp_b REAL(sp) :: pr, perc @@ -13049,8 +13053,12 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & INTEGER :: branch inv_cp = 1._sp/cp ps = cp*(1._sp-hp*hp)*TANH(pn*inv_cp)/(1._sp+hp*TANH(pn*inv_cp)) + CALL PUSHREAL4(ps) + ps = (1._sp+fq_ps)*ps es = hp*cp*(2._sp-hp)*TANH(en*inv_cp)/(1._sp+(1._sp-hp)*TANH(en*& & inv_cp)) + CALL PUSHREAL4(es) + es = (1._sp+fq_es)*es hp_imd = hp + (ps-es)*inv_cp IF (pn .GT. 0) THEN CALL PUSHCONTROL1B(0) @@ -13059,12 +13067,14 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & END IF pwx1 = 1._sp + (hp_imd/beta)**4 pwr1 = pwx1**(-0.25_sp) + CALL PUSHREAL4(perc) perc = hp_imd*cp*(1._sp-pwr1) pwx1 = 1._sp + (hp_imd/beta)**4 pwr1 = pwx1**(-0.25_sp) inv_cp = 1._sp/cp perc_b = perc_b - inv_cp*hp_b inv_cp_b = -(perc*hp_b) + CALL POPREAL4(perc) cp_b = cp_b + hp_imd*(1._sp-pwr1)*perc_b pwr1_b = -(hp_imd*cp*perc_b) pwx1_b = -(0.25_sp*pwx1**(-1.25)*pwr1_b) @@ -13072,15 +13082,18 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & & 4 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN - pn_b = pr_b + pn_b = pn_b + pr_b hp_imd_b = hp_imd_b - cp*pr_b hp_b = cp*pr_b cp_b = cp_b - (hp_imd-hp)*pr_b ELSE hp_b = 0.0_4 - pn_b = 0.0_4 END IF es_b = -(inv_cp*hp_imd_b) + inv_cp_b = inv_cp_b + (ps-es)*hp_imd_b + CALL POPREAL4(es) + fq_es_b = fq_es_b + es*es_b + es_b = (fq_es+1._sp)*es_b temp4 = TANH(en*inv_cp) temp3 = (-hp+1._sp)*temp4 + 1._sp temp1 = TANH(en*inv_cp) @@ -13093,8 +13106,11 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & ps_b = inv_cp*hp_imd_b temp_b4 = (1.0-TANH(en*inv_cp)**2)*temp0*temp_b3 temp_b5 = (1.0-TANH(en*inv_cp)**2)*(1._sp-hp)*temp_b0 - en_b = inv_cp*temp_b5 + inv_cp*temp_b4 + en_b = en_b + inv_cp*temp_b5 + inv_cp*temp_b4 cp_b = cp_b + hp*temp_b + CALL POPREAL4(ps) + fq_ps_b = fq_ps_b + ps*ps_b + ps_b = (fq_ps+1._sp)*ps_b temp = TANH(pn*inv_cp) temp0 = hp*temp + 1._sp temp1 = TANH(pn*inv_cp) @@ -13104,15 +13120,15 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & temp_b1 = -(temp2*temp1*temp_b/temp0) hp_b = hp_b + temp*temp_b1 - 2*hp*cp*temp1*temp_b temp_b2 = (1.0-TANH(pn*inv_cp)**2)*hp*temp_b1 - inv_cp_b = inv_cp_b + (ps-es)*hp_imd_b + en*temp_b5 + en*temp_b4 + & -& pn*temp_b2 + pn*temp_b0 + inv_cp_b = inv_cp_b + en*temp_b5 + en*temp_b4 + pn*temp_b2 + pn*& +& temp_b0 cp_b = cp_b + (1._sp-hp**2)*temp1*temp_b - inv_cp_b/cp**2 pn_b = pn_b + inv_cp*temp_b2 + inv_cp*temp_b0 END SUBROUTINE GR_PRODUCTION_B - SUBROUTINE GR_PRODUCTION(pn, en, cp, beta, hp, pr, perc) + SUBROUTINE GR_PRODUCTION(fq_ps, fq_es, pn, en, cp, beta, hp, pr, perc) IMPLICIT NONE - REAL(sp), INTENT(IN) :: pn, en, cp, beta + REAL(sp), INTENT(IN) :: fq_ps, fq_es, pn, en, cp, beta REAL(sp), INTENT(INOUT) :: hp REAL(sp), INTENT(OUT) :: pr, perc REAL(sp) :: inv_cp, ps, es, hp_imd @@ -13122,8 +13138,10 @@ SUBROUTINE GR_PRODUCTION(pn, en, cp, beta, hp, pr, perc) inv_cp = 1._sp/cp pr = 0._sp ps = cp*(1._sp-hp*hp)*TANH(pn*inv_cp)/(1._sp+hp*TANH(pn*inv_cp)) + ps = (1._sp+fq_ps)*ps es = hp*cp*(2._sp-hp)*TANH(en*inv_cp)/(1._sp+(1._sp-hp)*TANH(en*& & inv_cp)) + es = (1._sp+fq_es)*es hp_imd = hp + (ps-es)*inv_cp IF (pn .GT. 0) pr = pn - (hp_imd-hp)*cp pwx1 = 1._sp + (hp_imd/beta)**4 @@ -13134,42 +13152,46 @@ END SUBROUTINE GR_PRODUCTION ! Differentiation of gr_exchange in forward (tangent) mode (with options fixinterface noISIZE context): ! variations of useful results: l -! with respect to varying inputs: kexc ht - SUBROUTINE GR_EXCHANGE_D(kexc, kexc_d, ht, ht_d, l, l_d) +! with respect to varying inputs: kexc fq_l ht + SUBROUTINE GR_EXCHANGE_D(fq_l, fq_l_d, kexc, kexc_d, ht, ht_d, l, l_d) IMPLICIT NONE - REAL(sp), INTENT(IN) :: kexc - REAL(sp), INTENT(IN) :: kexc_d + REAL(sp), INTENT(IN) :: fq_l, kexc + REAL(sp), INTENT(IN) :: fq_l_d, kexc_d REAL(sp), INTENT(INOUT) :: ht REAL(sp), INTENT(INOUT) :: ht_d REAL(sp), INTENT(OUT) :: l REAL(sp), INTENT(OUT) :: l_d REAL(sp) :: temp temp = ht**3.5_sp - l_d = temp*kexc_d + kexc*3.5_sp*ht**2.5*ht_d - l = kexc*temp + l_d = temp*(kexc*fq_l_d+(fq_l+1._sp)*kexc_d) + (fq_l+1._sp)*kexc*& +& 3.5_sp*ht**2.5*ht_d + l = (fq_l+1._sp)*kexc*temp END SUBROUTINE GR_EXCHANGE_D ! Differentiation of gr_exchange in reverse (adjoint) mode (with options fixinterface noISIZE context): -! gradient of useful results: l kexc ht -! with respect to varying inputs: kexc ht - SUBROUTINE GR_EXCHANGE_B(kexc, kexc_b, ht, ht_b, l, l_b) +! gradient of useful results: l kexc fq_l ht +! with respect to varying inputs: kexc fq_l ht + SUBROUTINE GR_EXCHANGE_B(fq_l, fq_l_b, kexc, kexc_b, ht, ht_b, l, l_b) IMPLICIT NONE - REAL(sp), INTENT(IN) :: kexc - REAL(sp) :: kexc_b + REAL(sp), INTENT(IN) :: fq_l, kexc + REAL(sp) :: fq_l_b, kexc_b REAL(sp), INTENT(INOUT) :: ht REAL(sp), INTENT(INOUT) :: ht_b REAL(sp) :: l REAL(sp) :: l_b - kexc_b = kexc_b + ht**3.5_sp*l_b - ht_b = ht_b + 3.5_sp*ht**2.5*kexc*l_b + REAL(sp) :: temp_b + temp_b = ht**3.5_sp*l_b + ht_b = ht_b + 3.5_sp*ht**2.5*(fq_l+1._sp)*kexc*l_b + fq_l_b = fq_l_b + kexc*temp_b + kexc_b = kexc_b + (fq_l+1._sp)*temp_b END SUBROUTINE GR_EXCHANGE_B - SUBROUTINE GR_EXCHANGE(kexc, ht, l) + SUBROUTINE GR_EXCHANGE(fq_l, kexc, ht, l) IMPLICIT NONE - REAL(sp), INTENT(IN) :: kexc + REAL(sp), INTENT(IN) :: fq_l, kexc REAL(sp), INTENT(INOUT) :: ht REAL(sp), INTENT(OUT) :: l - l = kexc*ht**3.5_sp + l = (1._sp+fq_l)*kexc*ht**3.5_sp END SUBROUTINE GR_EXCHANGE ! Differentiation of gr_threshold_exchange in forward (tangent) mode (with options fixinterface noISIZE context): @@ -13553,610 +13575,78 @@ END SUBROUTINE GR_EXPONENTIAL_TRANSFER_D ! gradient of useful results: qe he be ! with respect to varying inputs: he be pre SUBROUTINE GR_EXPONENTIAL_TRANSFER_B(pre, pre_b, be, be_b, he, he_b, & -& qe, qe_b) - IMPLICIT NONE - REAL(sp), INTENT(IN) :: pre, be - REAL(sp) :: pre_b, be_b - REAL(sp), INTENT(INOUT) :: he - REAL(sp), INTENT(INOUT) :: he_b - REAL(sp) :: qe - REAL(sp) :: qe_b - REAL(sp) :: he_star, ar - REAL(sp) :: he_star_b, ar_b - INTRINSIC EXP - INTRINSIC LOG - REAL(sp) :: arg1 - REAL(sp) :: arg1_b - REAL(sp) :: temp - INTEGER :: branch - he_star = he + pre - ar = he_star/be - IF (ar .LT. -7._sp) THEN - CALL PUSHCONTROL2B(0) - ELSE IF (ar .GT. 7._sp) THEN - CALL PUSHCONTROL2B(1) - ELSE - arg1 = EXP(ar) + 1._sp - CALL PUSHCONTROL2B(2) - END IF - he_star_b = he_b - qe_b = qe_b - he_b - CALL POPCONTROL2B(branch) - IF (branch .EQ. 0) THEN - ar = he_star/be - be_b = be_b + EXP(ar)*qe_b - ar_b = EXP(ar)*be*qe_b - ELSE IF (branch .EQ. 1) THEN - ar = he_star/be - temp = EXP(ar) - he_star_b = he_star_b + qe_b - be_b = be_b + qe_b/temp - ar_b = -(EXP(ar)*be*qe_b/temp**2) - ELSE - be_b = be_b + LOG(arg1)*qe_b - arg1_b = be*qe_b/arg1 - ar = he_star/be - ar_b = EXP(ar)*arg1_b - END IF - he_star_b = he_star_b + ar_b/be - be_b = be_b - he_star*ar_b/be**2 - he_b = he_star_b - pre_b = he_star_b - END SUBROUTINE GR_EXPONENTIAL_TRANSFER_B - - SUBROUTINE GR_EXPONENTIAL_TRANSFER(pre, be, he, qe) - IMPLICIT NONE - REAL(sp), INTENT(IN) :: pre, be - REAL(sp), INTENT(INOUT) :: he - REAL(sp), INTENT(OUT) :: qe - REAL(sp) :: he_star, ar - INTRINSIC EXP - INTRINSIC LOG - REAL(sp) :: arg1 - he_star = he + pre - ar = he_star/be - IF (ar .LT. -7._sp) THEN - qe = be*EXP(ar) - ELSE IF (ar .GT. 7._sp) THEN - qe = he_star + be/EXP(ar) - ELSE - arg1 = EXP(ar) + 1._sp - qe = be*LOG(arg1) - END IF - he = he_star - qe - END SUBROUTINE GR_EXPONENTIAL_TRANSFER - -! Differentiation of gr_production_transfer_mlp_alg in forward (tangent) mode (with options fixinterface noISIZE context): -! variations of useful results: q hp ht -! with respect to varying inputs: kexc hp ht en f_q cp pn ct - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG_D(f_q, f_q_d, pn, pn_d, en, & -& en_d, cp, cp_d, beta, ct, ct_d, kexc, kexc_d, n, prcp, hp, hp_d, ht& -& , ht_d, q, q_d, pr, perc, l, prr, prd, qr, qd) - IMPLICIT NONE -! fixed NN output size - REAL(sp), DIMENSION(4), INTENT(IN) :: f_q - REAL(sp), DIMENSION(4), INTENT(IN) :: f_q_d - REAL(sp), INTENT(IN) :: pn, en, cp, beta, ct, kexc, n, prcp - REAL(sp), INTENT(IN) :: pn_d, en_d, cp_d, ct_d, kexc_d - REAL(sp), INTENT(INOUT) :: hp, ht, q - REAL(sp), INTENT(INOUT) :: hp_d, ht_d, q_d - REAL(sp), INTENT(OUT) :: pr, perc, l, prr, prd, qr, qd - REAL(sp) :: pr_d - REAL(sp) :: inv_cp, ps, es, hp_imd, pr_imd, ht_imd, nm1, d1pnm1 - REAL(sp) :: inv_cp_d, ps_d, es_d, hp_imd_d, pr_imd_d, ht_imd_d - INTRINSIC TANH - INTRINSIC MAX - REAL(sp) :: pwx1 - REAL(sp) :: pwx1_d - REAL(sp) :: pwr1 - REAL(sp) :: pwr1_d - REAL(sp) :: pwy1 - REAL(sp) :: pwy2 - REAL(sp) :: pwr2 - REAL(sp) :: pwr2_d - REAL(sp) :: pwx3 - REAL(sp) :: pwx3_d - REAL(sp) :: pwy3 - REAL(sp) :: pwr3 - REAL(sp) :: pwr3_d - REAL(sp) :: temp - REAL(sp) :: temp0 - REAL(sp) :: temp1 - REAL(sp) :: temp2 - REAL(sp) :: perc_d - REAL(sp) :: qr_d - REAL(sp) :: prd_d - REAL(sp) :: l_d - REAL(sp) :: qd_d - REAL(sp) :: prr_d -!% Production - inv_cp_d = -(cp_d/cp**2) - inv_cp = 1._sp/cp - pr = 0._sp - temp = TANH(pn*inv_cp) - temp0 = TANH(pn*inv_cp) - temp1 = cp*(-(hp*hp)+1._sp) - temp2 = temp1*temp0/(hp*temp+1._sp) - ps_d = (temp0*((1._sp-hp**2)*cp_d-cp*2*hp*hp_d)+temp1*(1.0-TANH(pn*& -& inv_cp)**2)*(inv_cp*pn_d+pn*inv_cp_d)-temp2*(temp*hp_d+hp*(1.0-& -& TANH(pn*inv_cp)**2)*(inv_cp*pn_d+pn*inv_cp_d)))/(hp*temp+1._sp) - ps = temp2 - ps_d = ps*f_q_d(1) + (f_q(1)+1._sp)*ps_d - ps = (1._sp+f_q(1))*ps - temp2 = TANH(en*inv_cp) - temp1 = TANH(en*inv_cp) - temp0 = hp*cp*(-hp+2._sp) - temp = temp0*temp1/((-hp+1._sp)*temp2+1._sp) - es_d = (temp1*((2._sp-hp)*(cp*hp_d+hp*cp_d)-hp*cp*hp_d)+temp0*(1.0-& -& TANH(en*inv_cp)**2)*(inv_cp*en_d+en*inv_cp_d)-temp*((1._sp-hp)*(& -& 1.0-TANH(en*inv_cp)**2)*(inv_cp*en_d+en*inv_cp_d)-temp2*hp_d))/((& -& 1._sp-hp)*temp2+1._sp) - es = temp - es_d = es*f_q_d(2) + (f_q(2)+1._sp)*es_d - es = (1._sp+f_q(2))*es - hp_imd_d = hp_d + inv_cp*(ps_d-es_d) + (ps-es)*inv_cp_d - hp_imd = hp + (ps-es)*inv_cp - IF (pn .GT. 0) THEN - pr_d = pn_d - cp*(hp_imd_d-hp_d) - (hp_imd-hp)*cp_d - pr = pn - (hp_imd-hp)*cp - ELSE - pr_d = 0.0_4 - END IF - pwx1_d = 4*hp_imd**3*hp_imd_d/beta**4 - pwx1 = 1._sp + (hp_imd/beta)**4 - pwr1_d = -(0.25_sp*pwx1**(-1.25)*pwx1_d) - pwr1 = pwx1**(-0.25_sp) - perc_d = (1._sp-pwr1)*(cp*hp_imd_d+hp_imd*cp_d) - hp_imd*cp*pwr1_d - perc = hp_imd*cp*(1._sp-pwr1) - hp_d = hp_imd_d - inv_cp*perc_d - perc*inv_cp_d - hp = hp_imd - perc*inv_cp -!% Transfer - nm1 = n - 1._sp - d1pnm1 = 1._sp/nm1 - temp2 = ht**3.5_sp - l_d = temp2*(kexc*f_q_d(4)+(f_q(4)+1._sp)*kexc_d) + (f_q(4)+1._sp)*& -& kexc*3.5_sp*ht**2.5*ht_d - l = (f_q(4)+1._sp)*kexc*temp2 - prr_d = 0.9_sp*((1._sp-f_q(3)**2)*(pr_d+perc_d)-(pr+perc)*2*f_q(3)*& -& f_q_d(3)) + l_d - prr = 0.9_sp*(1._sp-f_q(3)**2)*(pr+perc) + l - temp2 = 0.9_sp*(f_q(3)*f_q(3)) + 0.1_sp - prd_d = (pr+perc)*0.9_sp*2*f_q(3)*f_q_d(3) + temp2*(pr_d+perc_d) - prd = temp2*(pr+perc) - IF (prcp .LT. 0._sp) THEN - pwx1_d = ct*ht_d + ht*ct_d - pwx1 = ht*ct - pwy1 = -nm1 - IF (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1))) & -& THEN - pwr1_d = 0.0_4 - ELSE - pwr1_d = pwy1*pwx1**(pwy1-1)*pwx1_d - END IF - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - IF (ct .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2))) & -& THEN - pwr2_d = 0.0_4 - ELSE - pwr2_d = pwy2*ct**(pwy2-1)*ct_d - END IF - pwr2 = ct**pwy2 - pwx3_d = pwr1_d - pwr2_d - pwx3 = pwr1 - pwr2 - pwy3 = -d1pnm1 - IF (pwx3 .LE. 0.0 .AND. (pwy3 .EQ. 0.0 .OR. pwy3 .NE. INT(pwy3))) & -& THEN - pwr3_d = 0.0_4 - ELSE - pwr3_d = pwy3*pwx3**(pwy3-1)*pwx3_d - END IF - pwr3 = pwx3**pwy3 - pr_imd_d = pwr3_d - ct*ht_d - ht*ct_d - pr_imd = pwr3 - ht*ct - ELSE - pr_imd_d = prr_d - pr_imd = prr - END IF - IF (1.e-6_sp .LT. ht + pr_imd/ct) THEN - ht_imd_d = ht_d + (pr_imd_d-pr_imd*ct_d/ct)/ct - ht_imd = ht + pr_imd/ct - ELSE - ht_imd = 1.e-6_sp - ht_imd_d = 0.0_4 - END IF - pwx1_d = ct*ht_imd_d + ht_imd*ct_d - pwx1 = ht_imd*ct - pwy1 = -nm1 - IF (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1))) & -& THEN - pwr1_d = 0.0_4 - ELSE - pwr1_d = pwy1*pwx1**(pwy1-1)*pwx1_d - END IF - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - IF (ct .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2))) THEN - pwr2_d = 0.0_4 - ELSE - pwr2_d = pwy2*ct**(pwy2-1)*ct_d - END IF - pwr2 = ct**pwy2 - pwx3_d = pwr1_d + pwr2_d - pwx3 = pwr1 + pwr2 - pwy3 = -d1pnm1 - IF (pwx3 .LE. 0.0 .AND. (pwy3 .EQ. 0.0 .OR. pwy3 .NE. INT(pwy3))) & -& THEN - pwr3_d = 0.0_4 - ELSE - pwr3_d = pwy3*pwx3**(pwy3-1)*pwx3_d - END IF - pwr3 = pwx3**pwy3 - ht_d = (pwr3_d-pwr3*ct_d/ct)/ct - ht = pwr3/ct - IF (0._sp .LT. prd + l) THEN - qd_d = prd_d + l_d - qd = prd + l - ELSE - qd = 0._sp - qd_d = 0.0_4 - END IF - qr_d = ct*(ht_imd_d-ht_d) + (ht_imd-ht)*ct_d - qr = (ht_imd-ht)*ct - q_d = qr_d + qd_d - q = qr + qd - END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG_D - -! Differentiation of gr_production_transfer_mlp_alg in reverse (adjoint) mode (with options fixinterface noISIZE context): -! gradient of useful results: q kexc hp ht en f_q cp pn ct -! with respect to varying inputs: kexc hp ht en f_q cp pn ct - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG_B(f_q, f_q_b, pn, pn_b, en, & -& en_b, cp, cp_b, beta, ct, ct_b, kexc, kexc_b, n, prcp, hp, hp_b, ht& -& , ht_b, q, q_b, pr, perc, l, prr, prd, qr, qd) - IMPLICIT NONE -! fixed NN output size - REAL(sp), DIMENSION(4), INTENT(IN) :: f_q - REAL(sp), DIMENSION(4) :: f_q_b - REAL(sp), INTENT(IN) :: pn, en, cp, beta, ct, kexc, n, prcp - REAL(sp) :: pn_b, en_b, cp_b, ct_b, kexc_b - REAL(sp), INTENT(INOUT) :: hp, ht, q - REAL(sp), INTENT(INOUT) :: hp_b, ht_b, q_b - REAL(sp) :: pr, perc, l, prr, prd, qr, qd - REAL(sp) :: inv_cp, ps, es, hp_imd, pr_imd, ht_imd, nm1, d1pnm1 - REAL(sp) :: inv_cp_b, ps_b, es_b, hp_imd_b, pr_imd_b, ht_imd_b - INTRINSIC TANH - INTRINSIC MAX - REAL(sp) :: pwx1 - REAL(sp) :: pwx1_b - REAL(sp) :: pwr1 - REAL(sp) :: pwr1_b - REAL(sp) :: pwy1 - REAL(sp) :: pwy2 - REAL(sp) :: pwr2 - REAL(sp) :: pwr2_b - REAL(sp) :: pwx3 - REAL(sp) :: pwx3_b - REAL(sp) :: pwy3 - REAL(sp) :: pwr3 - REAL(sp) :: pwr3_b - REAL(sp) :: temp - REAL(sp) :: temp0 - REAL(sp) :: temp_b - REAL(sp) :: temp1 - REAL(sp) :: temp2 - REAL(sp) :: temp3 - REAL(sp) :: temp_b0 - REAL(sp) :: temp_b1 - REAL(sp) :: temp4 - REAL(sp) :: temp_b2 - REAL(sp) :: temp_b3 - REAL(sp) :: temp_b4 - REAL(sp) :: temp_b5 - INTEGER :: branch - REAL(sp) :: perc_b - REAL(sp) :: pr_b - REAL(sp) :: qr_b - REAL(sp) :: prd_b - REAL(sp) :: l_b - REAL(sp) :: qd_b - REAL(sp) :: prr_b -!% Production - inv_cp = 1._sp/cp - pr = 0._sp - ps = cp*(1._sp-hp*hp)*TANH(pn*inv_cp)/(1._sp+hp*TANH(pn*inv_cp)) - CALL PUSHREAL4(ps) - ps = (1._sp+f_q(1))*ps - es = hp*cp*(2._sp-hp)*TANH(en*inv_cp)/(1._sp+(1._sp-hp)*TANH(en*& -& inv_cp)) - CALL PUSHREAL4(es) - es = (1._sp+f_q(2))*es - hp_imd = hp + (ps-es)*inv_cp - IF (pn .GT. 0) THEN - pr = pn - (hp_imd-hp)*cp - CALL PUSHCONTROL1B(0) - ELSE - CALL PUSHCONTROL1B(1) - END IF - pwx1 = 1._sp + (hp_imd/beta)**4 - pwr1 = pwx1**(-0.25_sp) - perc = hp_imd*cp*(1._sp-pwr1) -!% Transfer - nm1 = n - 1._sp - d1pnm1 = 1._sp/nm1 - l = (1._sp+f_q(4))*kexc*ht**3.5_sp - prr = 0.9_sp*(1._sp-f_q(3)**2)*(pr+perc) + l - prd = (0.1_sp+0.9_sp*f_q(3)**2)*(pr+perc) - IF (prcp .LT. 0._sp) THEN - pwx1 = ht*ct - pwy1 = -nm1 - CALL PUSHREAL4(pwr1) - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - pwr2 = ct**pwy2 - pwx3 = pwr1 - pwr2 - pwy3 = -d1pnm1 - pwr3 = pwx3**pwy3 - pr_imd = pwr3 - ht*ct - CALL PUSHCONTROL1B(1) - ELSE - pr_imd = prr - CALL PUSHCONTROL1B(0) - END IF - IF (1.e-6_sp .LT. ht + pr_imd/ct) THEN - ht_imd = ht + pr_imd/ct - CALL PUSHCONTROL1B(0) - ELSE - ht_imd = 1.e-6_sp - CALL PUSHCONTROL1B(1) - END IF - CALL PUSHREAL4(pwx1) - pwx1 = ht_imd*ct - CALL PUSHREAL4(pwy1) - pwy1 = -nm1 - CALL PUSHREAL4(pwr1) - pwr1 = pwx1**pwy1 - CALL PUSHREAL4(pwy2) - pwy2 = -nm1 - pwr2 = ct**pwy2 - CALL PUSHREAL4(pwx3) - pwx3 = pwr1 + pwr2 - CALL PUSHREAL4(pwy3) - pwy3 = -d1pnm1 - pwr3 = pwx3**pwy3 - CALL PUSHREAL4(ht) - ht = pwr3/ct - IF (0._sp .LT. prd + l) THEN - CALL PUSHCONTROL1B(0) - ELSE - CALL PUSHCONTROL1B(1) - END IF - pwx1 = ht_imd*ct - nm1 = n - 1._sp - pwy1 = -nm1 - pwy2 = -nm1 - d1pnm1 = 1._sp/nm1 - pwy3 = -d1pnm1 - inv_cp = 1._sp/cp - qr_b = q_b - qd_b = q_b - ht_imd_b = ct*qr_b - ht_b = ht_b - ct*qr_b - ct_b = ct_b + (ht_imd-ht)*qr_b - CALL POPCONTROL1B(branch) - IF (branch .EQ. 0) THEN - prd_b = qd_b - l_b = qd_b - ELSE - l_b = 0.0_4 - prd_b = 0.0_4 - END IF - CALL POPREAL4(ht) - pwr3_b = ht_b/ct - IF (pwx3 .LE. 0.0 .AND. (pwy3 .EQ. 0.0 .OR. pwy3 .NE. INT(pwy3))) & -& THEN - pwx3_b = 0.0_4 - ELSE - pwx3_b = pwy3*pwx3**(pwy3-1)*pwr3_b - END IF - CALL POPREAL4(pwy3) - CALL POPREAL4(pwx3) - pwr1_b = pwx3_b - pwr2_b = pwx3_b - IF (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1))) & -& THEN - pwx1_b = 0.0_4 - ELSE - pwx1_b = pwy1*pwx1**(pwy1-1)*pwr1_b - END IF - IF (ct .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2))) THEN - ct_b = ct_b + ht_imd*pwx1_b - pwr3*ht_b/ct**2 - ELSE - ct_b = ct_b + pwy2*ct**(pwy2-1)*pwr2_b - pwr3*ht_b/ct**2 + ht_imd*& -& pwx1_b - END IF - CALL POPREAL4(pwy2) - CALL POPREAL4(pwr1) - CALL POPREAL4(pwy1) - CALL POPREAL4(pwx1) - ht_imd_b = ht_imd_b + ct*pwx1_b - CALL POPCONTROL1B(branch) - IF (branch .EQ. 0) THEN - ht_b = ht_imd_b - pr_imd_b = ht_imd_b/ct - ct_b = ct_b - pr_imd*ht_imd_b/ct**2 - ELSE - ht_b = 0.0_4 - pr_imd_b = 0.0_4 - END IF - CALL POPCONTROL1B(branch) - IF (branch .EQ. 0) THEN - prr_b = pr_imd_b - ELSE - pwr3_b = pr_imd_b - IF (pwx3 .LE. 0.0 .AND. (pwy3 .EQ. 0.0 .OR. pwy3 .NE. INT(pwy3))) & -& THEN - pwx3_b = 0.0_4 - ELSE - pwx3_b = pwy3*pwx3**(pwy3-1)*pwr3_b - END IF - pwr1_b = pwx3_b - pwr2_b = -pwx3_b - CALL POPREAL4(pwr1) - IF (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1))) & -& THEN - pwx1_b = 0.0_4 - ELSE - pwx1_b = pwy1*pwx1**(pwy1-1)*pwr1_b - END IF - ht_b = ht_b + ct*pwx1_b - ct*pr_imd_b - IF (ct .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2))) & -& THEN - ct_b = ct_b + ht*pwx1_b - ht*pr_imd_b - ELSE - ct_b = ct_b + pwy2*ct**(pwy2-1)*pwr2_b - ht*pr_imd_b + ht*pwx1_b - END IF - pwx1 = 1._sp + (hp_imd/beta)**4 - prr_b = 0.0_4 - END IF - f_q_b(3) = f_q_b(3) + 2*f_q(3)*0.9_sp*(pr+perc)*prd_b - 2*f_q(3)*(pr& -& +perc)*0.9_sp*prr_b - temp_b5 = (0.9_sp*f_q(3)**2+0.1_sp)*prd_b - pr_b = temp_b5 - perc_b = temp_b5 - temp_b5 = (1._sp-f_q(3)**2)*0.9_sp*prr_b - l_b = l_b + prr_b - pr_b = pr_b + temp_b5 - perc_b = perc_b + temp_b5 - inv_cp*hp_b - temp_b5 = ht**3.5_sp*l_b - ht_b = ht_b + 3.5_sp*ht**2.5*(f_q(4)+1._sp)*kexc*l_b - f_q_b(4) = f_q_b(4) + kexc*temp_b5 - kexc_b = kexc_b + (f_q(4)+1._sp)*temp_b5 - inv_cp_b = -(perc*hp_b) - cp_b = cp_b + hp_imd*(1._sp-pwr1)*perc_b - pwr1_b = -(hp_imd*cp*perc_b) - pwx1_b = -(0.25_sp*pwx1**(-1.25)*pwr1_b) - hp_imd_b = hp_b + cp*(1._sp-pwr1)*perc_b + 4*hp_imd**3*pwx1_b/beta**& -& 4 - CALL POPCONTROL1B(branch) - IF (branch .EQ. 0) THEN - pn_b = pn_b + pr_b - hp_imd_b = hp_imd_b - cp*pr_b - hp_b = cp*pr_b - cp_b = cp_b - (hp_imd-hp)*pr_b - ELSE - hp_b = 0.0_4 - END IF - es_b = -(inv_cp*hp_imd_b) - inv_cp_b = inv_cp_b + (ps-es)*hp_imd_b - CALL POPREAL4(es) - f_q_b(2) = f_q_b(2) + es*es_b - es_b = (f_q(2)+1._sp)*es_b - temp4 = TANH(en*inv_cp) - temp3 = (-hp+1._sp)*temp4 + 1._sp - temp1 = TANH(en*inv_cp) - temp0 = hp*cp*(-hp+2._sp) - temp_b3 = es_b/temp3 - temp_b = (2._sp-hp)*temp1*temp_b3 - temp_b0 = -(temp0*temp1*temp_b3/temp3) - hp_b = hp_b + hp_imd_b + cp*temp_b - hp*cp*temp1*temp_b3 - temp4*& -& temp_b0 - ps_b = inv_cp*hp_imd_b - temp_b4 = (1.0-TANH(en*inv_cp)**2)*temp0*temp_b3 - temp_b5 = (1.0-TANH(en*inv_cp)**2)*(1._sp-hp)*temp_b0 - en_b = en_b + inv_cp*temp_b5 + inv_cp*temp_b4 - cp_b = cp_b + hp*temp_b - CALL POPREAL4(ps) - f_q_b(1) = f_q_b(1) + ps*ps_b - ps_b = (f_q(1)+1._sp)*ps_b - temp = TANH(pn*inv_cp) - temp0 = hp*temp + 1._sp - temp1 = TANH(pn*inv_cp) - temp2 = cp*(-(hp*hp)+1._sp) - temp_b = ps_b/temp0 - temp_b0 = (1.0-TANH(pn*inv_cp)**2)*temp2*temp_b - temp_b1 = -(temp2*temp1*temp_b/temp0) - hp_b = hp_b + temp*temp_b1 - 2*hp*cp*temp1*temp_b - temp_b2 = (1.0-TANH(pn*inv_cp)**2)*hp*temp_b1 - inv_cp_b = inv_cp_b + en*temp_b5 + en*temp_b4 + pn*temp_b2 + pn*& -& temp_b0 - cp_b = cp_b + (1._sp-hp**2)*temp1*temp_b - inv_cp_b/cp**2 - pn_b = pn_b + inv_cp*temp_b2 + inv_cp*temp_b0 - END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG_B - - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG(f_q, pn, en, cp, beta, ct, & -& kexc, n, prcp, hp, ht, q, pr, perc, l, prr, prd, qr, qd) +& qe, qe_b) IMPLICIT NONE -! fixed NN output size - REAL(sp), DIMENSION(4), INTENT(IN) :: f_q - REAL(sp), INTENT(IN) :: pn, en, cp, beta, ct, kexc, n, prcp - REAL(sp), INTENT(INOUT) :: hp, ht, q - REAL(sp), INTENT(OUT) :: pr, perc, l, prr, prd, qr, qd - REAL(sp) :: inv_cp, ps, es, hp_imd, pr_imd, ht_imd, nm1, d1pnm1 - INTRINSIC TANH - INTRINSIC MAX - REAL(sp) :: pwx1 - REAL(sp) :: pwr1 - REAL(sp) :: pwy1 - REAL(sp) :: pwy2 - REAL(sp) :: pwr2 - REAL(sp) :: pwx3 - REAL(sp) :: pwy3 - REAL(sp) :: pwr3 -!% Production - inv_cp = 1._sp/cp - pr = 0._sp - ps = cp*(1._sp-hp*hp)*TANH(pn*inv_cp)/(1._sp+hp*TANH(pn*inv_cp)) - ps = (1._sp+f_q(1))*ps - es = hp*cp*(2._sp-hp)*TANH(en*inv_cp)/(1._sp+(1._sp-hp)*TANH(en*& -& inv_cp)) - es = (1._sp+f_q(2))*es - hp_imd = hp + (ps-es)*inv_cp - IF (pn .GT. 0) pr = pn - (hp_imd-hp)*cp - pwx1 = 1._sp + (hp_imd/beta)**4 - pwr1 = pwx1**(-0.25_sp) - perc = hp_imd*cp*(1._sp-pwr1) - hp = hp_imd - perc*inv_cp -!% Transfer - nm1 = n - 1._sp - d1pnm1 = 1._sp/nm1 - l = (1._sp+f_q(4))*kexc*ht**3.5_sp - prr = 0.9_sp*(1._sp-f_q(3)**2)*(pr+perc) + l - prd = (0.1_sp+0.9_sp*f_q(3)**2)*(pr+perc) - IF (prcp .LT. 0._sp) THEN - pwx1 = ht*ct - pwy1 = -nm1 - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - pwr2 = ct**pwy2 - pwx3 = pwr1 - pwr2 - pwy3 = -d1pnm1 - pwr3 = pwx3**pwy3 - pr_imd = pwr3 - ht*ct + REAL(sp), INTENT(IN) :: pre, be + REAL(sp) :: pre_b, be_b + REAL(sp), INTENT(INOUT) :: he + REAL(sp), INTENT(INOUT) :: he_b + REAL(sp) :: qe + REAL(sp) :: qe_b + REAL(sp) :: he_star, ar + REAL(sp) :: he_star_b, ar_b + INTRINSIC EXP + INTRINSIC LOG + REAL(sp) :: arg1 + REAL(sp) :: arg1_b + REAL(sp) :: temp + INTEGER :: branch + he_star = he + pre + ar = he_star/be + IF (ar .LT. -7._sp) THEN + CALL PUSHCONTROL2B(0) + ELSE IF (ar .GT. 7._sp) THEN + CALL PUSHCONTROL2B(1) ELSE - pr_imd = prr + arg1 = EXP(ar) + 1._sp + CALL PUSHCONTROL2B(2) END IF - IF (1.e-6_sp .LT. ht + pr_imd/ct) THEN - ht_imd = ht + pr_imd/ct + he_star_b = he_b + qe_b = qe_b - he_b + CALL POPCONTROL2B(branch) + IF (branch .EQ. 0) THEN + ar = he_star/be + be_b = be_b + EXP(ar)*qe_b + ar_b = EXP(ar)*be*qe_b + ELSE IF (branch .EQ. 1) THEN + ar = he_star/be + temp = EXP(ar) + he_star_b = he_star_b + qe_b + be_b = be_b + qe_b/temp + ar_b = -(EXP(ar)*be*qe_b/temp**2) ELSE - ht_imd = 1.e-6_sp + be_b = be_b + LOG(arg1)*qe_b + arg1_b = be*qe_b/arg1 + ar = he_star/be + ar_b = EXP(ar)*arg1_b END IF - pwx1 = ht_imd*ct - pwy1 = -nm1 - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - pwr2 = ct**pwy2 - pwx3 = pwr1 + pwr2 - pwy3 = -d1pnm1 - pwr3 = pwx3**pwy3 - ht = pwr3/ct - IF (0._sp .LT. prd + l) THEN - qd = prd + l + he_star_b = he_star_b + ar_b/be + be_b = be_b - he_star*ar_b/be**2 + he_b = he_star_b + pre_b = he_star_b + END SUBROUTINE GR_EXPONENTIAL_TRANSFER_B + + SUBROUTINE GR_EXPONENTIAL_TRANSFER(pre, be, he, qe) + IMPLICIT NONE + REAL(sp), INTENT(IN) :: pre, be + REAL(sp), INTENT(INOUT) :: he + REAL(sp), INTENT(OUT) :: qe + REAL(sp) :: he_star, ar + INTRINSIC EXP + INTRINSIC LOG + REAL(sp) :: arg1 + he_star = he + pre + ar = he_star/be + IF (ar .LT. -7._sp) THEN + qe = be*EXP(ar) + ELSE IF (ar .GT. 7._sp) THEN + qe = he_star + be/EXP(ar) ELSE - qd = 0._sp + arg1 = EXP(ar) + 1._sp + qe = be*LOG(arg1) END IF - qr = (ht_imd-ht)*ct - q = qr + qd - END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG + he = he_star - qe + END SUBROUTINE GR_EXPONENTIAL_TRANSFER ! Differentiation of gr_production_transfer_ode in forward (tangent) mode (with options fixinterface noISIZE context): ! variations of useful results: q hp ht @@ -14488,14 +13978,14 @@ END SUBROUTINE GR_PRODUCTION_TRANSFER_ODE ! Differentiation of gr_production_transfer_mlp_ode in forward (tangent) mode (with options fixinterface noISIZE context): ! variations of useful results: q hp ht -! with respect to varying inputs: kexc hp ht en f_q cp pn ct - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(f_q, f_q_d, pn, pn_d, en, & +! with respect to varying inputs: kexc hp ht en fq cp pn ct + SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(fq, fq_d, pn, pn_d, en, & & en_d, cp, cp_d, ct, ct_d, kexc, kexc_d, hp, hp_d, ht, ht_d, q, q_d, & & l) IMPLICIT NONE ! fixed NN output size - REAL(sp), DIMENSION(5), INTENT(IN) :: f_q - REAL(sp), DIMENSION(5), INTENT(IN) :: f_q_d + REAL(sp), DIMENSION(5), INTENT(IN) :: fq + REAL(sp), DIMENSION(5), INTENT(IN) :: fq_d REAL(sp), INTENT(IN) :: pn, en, cp, ct, kexc REAL(sp), INTENT(IN) :: pn_d, en_d, cp_d, ct_d, kexc_d REAL(sp), INTENT(INOUT) :: hp, ht, q @@ -14514,11 +14004,11 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(f_q, f_q_d, pn, pn_d, en, & ! dt = 1._sp/real(n_subtimesteps, sp) dt = 1._sp !do i = 1, n_subtimesteps - temp = (f_q(2)+1._sp)*(-hp+2._sp) - fhp_d = (1._sp-hp**2)*(pn*f_q_d(1)+(f_q(1)+1._sp)*pn_d) - (f_q(1)+& -& 1._sp)*pn*2*hp*hp_d - en*hp*((2._sp-hp)*f_q_d(2)-(f_q(2)+1._sp)*& -& hp_d) - temp*(hp*en_d+en*hp_d) - fhp = (f_q(1)+1._sp)*pn*(1._sp-hp*hp) - temp*(en*hp) + temp = (fq(2)+1._sp)*(-hp+2._sp) + fhp_d = (1._sp-hp**2)*(pn*fq_d(1)+(fq(1)+1._sp)*pn_d) - (fq(1)+1._sp& +& )*pn*2*hp*hp_d - en*hp*((2._sp-hp)*fq_d(2)-(fq(2)+1._sp)*hp_d) - & +& temp*(hp*en_d+en*hp_d) + fhp = (fq(1)+1._sp)*pn*(1._sp-hp*hp) - temp*(en*hp) hp_d = hp_d + dt*(inv_cp*fhp_d+fhp*inv_cp_d) hp = hp + dt*fhp*inv_cp IF (hp .LE. 0._sp) THEN @@ -14529,16 +14019,16 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(f_q, f_q_d, pn, pn_d, en, & hp = 1._sp - 1.e-6_sp hp_d = 0.0_4 END IF - temp = (-(f_q(3)*f_q(3))+1._sp)*(hp*hp) + temp = (-(fq(3)*fq(3))+1._sp)*(hp*hp) temp0 = ht**3.5_sp temp1 = ht**5 - fht_d = 0.9_sp*((f_q(1)+1._sp)*pn*((1._sp-f_q(3)**2)*2*hp*hp_d-hp**2& -& *2*f_q(3)*f_q_d(3))+temp*(pn*f_q_d(1)+(f_q(1)+1._sp)*pn_d)) + & -& temp0*(kexc*f_q_d(4)+(f_q(4)+1._sp)*kexc_d) + ((f_q(4)+1._sp)*kexc& -& *3.5_sp*ht**2.5-(f_q(5)+1._sp)*ct*5*ht**4)*ht_d - temp1*(ct*f_q_d(& -& 5)+(f_q(5)+1._sp)*ct_d) - fht = 0.9_sp*(temp*((f_q(1)+1._sp)*pn)) + (f_q(4)+1._sp)*kexc*temp0 & -& - (f_q(5)+1._sp)*ct*temp1 + fht_d = 0.9_sp*((fq(1)+1._sp)*pn*((1._sp-fq(3)**2)*2*hp*hp_d-hp**2*2& +& *fq(3)*fq_d(3))+temp*(pn*fq_d(1)+(fq(1)+1._sp)*pn_d)) + temp0*(& +& kexc*fq_d(4)+(fq(4)+1._sp)*kexc_d) + ((fq(4)+1._sp)*kexc*3.5_sp*ht& +& **2.5-(fq(5)+1._sp)*ct*5*ht**4)*ht_d - temp1*(ct*fq_d(5)+(fq(5)+& +& 1._sp)*ct_d) + fht = 0.9_sp*(temp*((fq(1)+1._sp)*pn)) + (fq(4)+1._sp)*kexc*temp0 - & +& (fq(5)+1._sp)*ct*temp1 ht_d = ht_d + dt*(fht_d-fht*ct_d/ct)/ct ht = ht + dt*fht/ct IF (ht .LE. 0._sp) THEN @@ -14551,28 +14041,28 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(f_q, f_q_d, pn, pn_d, en, & END IF !end do temp1 = ht**3.5_sp - l_d = temp1*(kexc*f_q_d(4)+(f_q(4)+1._sp)*kexc_d) + (f_q(4)+1._sp)*& -& kexc*3.5_sp*ht**2.5*ht_d - l = (f_q(4)+1._sp)*kexc*temp1 + l_d = temp1*(kexc*fq_d(4)+(fq(4)+1._sp)*kexc_d) + (fq(4)+1._sp)*kexc& +& *3.5_sp*ht**2.5*ht_d + l = (fq(4)+1._sp)*kexc*temp1 temp1 = ht**5 - temp0 = (f_q(1)+1._sp)*pn*(hp*hp) - temp = 0.9_sp*(f_q(3)*f_q(3)) + 0.1_sp - q_d = temp1*(ct*f_q_d(5)+(f_q(5)+1._sp)*ct_d) + (f_q(5)+1._sp)*ct*5*& -& ht**4*ht_d + temp0*0.9_sp*2*f_q(3)*f_q_d(3) + temp*(hp**2*(pn*& -& f_q_d(1)+(f_q(1)+1._sp)*pn_d)+(f_q(1)+1._sp)*pn*2*hp*hp_d) + l_d - q = (f_q(5)+1._sp)*ct*temp1 + temp*temp0 + l + temp0 = (fq(1)+1._sp)*pn*(hp*hp) + temp = 0.9_sp*(fq(3)*fq(3)) + 0.1_sp + q_d = temp1*(ct*fq_d(5)+(fq(5)+1._sp)*ct_d) + (fq(5)+1._sp)*ct*5*ht& +& **4*ht_d + temp0*0.9_sp*2*fq(3)*fq_d(3) + temp*(hp**2*(pn*fq_d(1)+& +& (fq(1)+1._sp)*pn_d)+(fq(1)+1._sp)*pn*2*hp*hp_d) + l_d + q = (fq(5)+1._sp)*ct*temp1 + temp*temp0 + l END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D ! Differentiation of gr_production_transfer_mlp_ode in reverse (adjoint) mode (with options fixinterface noISIZE context): -! gradient of useful results: q kexc hp ht en f_q cp pn ct -! with respect to varying inputs: kexc hp ht en f_q cp pn ct - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & +! gradient of useful results: q kexc hp ht en fq cp pn ct +! with respect to varying inputs: kexc hp ht en fq cp pn ct + SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(fq, fq_b, pn, pn_b, en, & & en_b, cp, cp_b, ct, ct_b, kexc, kexc_b, hp, hp_b, ht, ht_b, q, q_b, & & l) IMPLICIT NONE ! fixed NN output size - REAL(sp), DIMENSION(5), INTENT(IN) :: f_q - REAL(sp), DIMENSION(5) :: f_q_b + REAL(sp), DIMENSION(5), INTENT(IN) :: fq + REAL(sp), DIMENSION(5) :: fq_b REAL(sp), INTENT(IN) :: pn, en, cp, ct, kexc REAL(sp) :: pn_b, en_b, cp_b, ct_b, kexc_b REAL(sp), INTENT(INOUT) :: hp, ht, q @@ -14592,8 +14082,8 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & ! dt = 1._sp/real(n_subtimesteps, sp) dt = 1._sp !do i = 1, n_subtimesteps - fhp = (1._sp+f_q(1))*pn*(1._sp-hp**2) - (1._sp+f_q(2))*en*hp*(2._sp-& -& hp) + fhp = (1._sp+fq(1))*pn*(1._sp-hp**2) - (1._sp+fq(2))*en*hp*(2._sp-hp& +& ) CALL PUSHREAL4(hp) hp = hp + dt*fhp*inv_cp IF (hp .LE. 0._sp) THEN @@ -14608,8 +14098,8 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & ELSE CALL PUSHCONTROL1B(1) END IF - fht = 0.9_sp*(1._sp-f_q(3)**2)*(1._sp+f_q(1))*pn*hp**2 + (1._sp+f_q(& -& 4))*kexc*ht**3.5_sp - (1._sp+f_q(5))*ct*ht**5 + fht = 0.9_sp*(1._sp-fq(3)**2)*(1._sp+fq(1))*pn*hp**2 + (1._sp+fq(4))& +& *kexc*ht**3.5_sp - (1._sp+fq(5))*ct*ht**5 CALL PUSHREAL4(ht) ht = ht + dt*fht/ct IF (ht .LE. 0._sp) THEN @@ -14626,19 +14116,19 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & END IF l_b = q_b temp_b2 = ht**5*q_b - ht_b = ht_b + 5*ht**4*(f_q(5)+1._sp)*ct*q_b + 3.5_sp*ht**2.5*(f_q(4)& -& +1._sp)*kexc*l_b - f_q_b(3) = f_q_b(3) + 2*f_q(3)*0.9_sp*(f_q(1)+1._sp)*pn*hp**2*q_b - temp_b1 = (0.9_sp*f_q(3)**2+0.1_sp)*q_b + ht_b = ht_b + 5*ht**4*(fq(5)+1._sp)*ct*q_b + 3.5_sp*ht**2.5*(fq(4)+& +& 1._sp)*kexc*l_b + fq_b(3) = fq_b(3) + 2*fq(3)*0.9_sp*(fq(1)+1._sp)*pn*hp**2*q_b + temp_b1 = (0.9_sp*fq(3)**2+0.1_sp)*q_b temp_b0 = hp**2*temp_b1 - hp_b = hp_b + 2*hp*(f_q(1)+1._sp)*pn*temp_b1 - f_q_b(1) = f_q_b(1) + pn*temp_b0 - pn_b = pn_b + (f_q(1)+1._sp)*temp_b0 - f_q_b(5) = f_q_b(5) + ct*temp_b2 - ct_b = ct_b + (f_q(5)+1._sp)*temp_b2 + hp_b = hp_b + 2*hp*(fq(1)+1._sp)*pn*temp_b1 + fq_b(1) = fq_b(1) + pn*temp_b0 + pn_b = pn_b + (fq(1)+1._sp)*temp_b0 + fq_b(5) = fq_b(5) + ct*temp_b2 + ct_b = ct_b + (fq(5)+1._sp)*temp_b2 temp_b2 = ht**3.5_sp*l_b - f_q_b(4) = f_q_b(4) + kexc*temp_b2 - kexc_b = kexc_b + (f_q(4)+1._sp)*temp_b2 + fq_b(4) = fq_b(4) + kexc*temp_b2 + kexc_b = kexc_b + (fq(4)+1._sp)*temp_b2 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) ht_b = 0.0_4 CALL POPCONTROL1B(branch) @@ -14648,20 +14138,20 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & temp_b2 = dt*ht_b/ct fht_b = temp_b2 ct_b = ct_b - fht*temp_b2/ct - temp_b1 = (f_q(1)+1._sp)*pn*0.9_sp*fht_b - temp_b0 = (1._sp-f_q(3)**2)*hp**2*0.9_sp*fht_b + temp_b1 = (fq(1)+1._sp)*pn*0.9_sp*fht_b + temp_b0 = (1._sp-fq(3)**2)*hp**2*0.9_sp*fht_b temp_b = ht**3.5_sp*fht_b - ht_b = ht_b + (3.5_sp*ht**2.5*(f_q(4)+1._sp)*kexc-5*ht**4*(f_q(5)+& + ht_b = ht_b + (3.5_sp*ht**2.5*(fq(4)+1._sp)*kexc-5*ht**4*(fq(5)+& & 1._sp)*ct)*fht_b temp_b2 = -(ht**5*fht_b) - f_q_b(5) = f_q_b(5) + ct*temp_b2 - ct_b = ct_b + (f_q(5)+1._sp)*temp_b2 - f_q_b(4) = f_q_b(4) + kexc*temp_b - kexc_b = kexc_b + (f_q(4)+1._sp)*temp_b - f_q_b(1) = f_q_b(1) + pn*temp_b0 - pn_b = pn_b + (f_q(1)+1._sp)*temp_b0 - f_q_b(3) = f_q_b(3) - 2*f_q(3)*hp**2*temp_b1 - hp_b = hp_b + 2*hp*(1._sp-f_q(3)**2)*temp_b1 + fq_b(5) = fq_b(5) + ct*temp_b2 + ct_b = ct_b + (fq(5)+1._sp)*temp_b2 + fq_b(4) = fq_b(4) + kexc*temp_b + kexc_b = kexc_b + (fq(4)+1._sp)*temp_b + fq_b(1) = fq_b(1) + pn*temp_b0 + pn_b = pn_b + (fq(1)+1._sp)*temp_b0 + fq_b(3) = fq_b(3) - 2*fq(3)*hp**2*temp_b1 + hp_b = hp_b + 2*hp*(1._sp-fq(3)**2)*temp_b1 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) hp_b = 0.0_4 CALL POPCONTROL1B(branch) @@ -14672,21 +14162,21 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & inv_cp_b = fhp*dt*hp_b temp_b = (1._sp-hp**2)*fhp_b temp_b0 = -(en*hp*fhp_b) - temp_b1 = -((f_q(2)+1._sp)*(2._sp-hp)*fhp_b) - hp_b = hp_b + en*temp_b1 - 2*hp*(f_q(1)+1._sp)*pn*fhp_b - (f_q(2)+& + temp_b1 = -((fq(2)+1._sp)*(2._sp-hp)*fhp_b) + hp_b = hp_b + en*temp_b1 - 2*hp*(fq(1)+1._sp)*pn*fhp_b - (fq(2)+& & 1._sp)*temp_b0 en_b = en_b + hp*temp_b1 - f_q_b(2) = f_q_b(2) + (2._sp-hp)*temp_b0 - f_q_b(1) = f_q_b(1) + pn*temp_b - pn_b = pn_b + (f_q(1)+1._sp)*temp_b + fq_b(2) = fq_b(2) + (2._sp-hp)*temp_b0 + fq_b(1) = fq_b(1) + pn*temp_b + pn_b = pn_b + (fq(1)+1._sp)*temp_b cp_b = cp_b - inv_cp_b/cp**2 END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE(f_q, pn, en, cp, ct, kexc, & -& hp, ht, q, l) + SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE(fq, pn, en, cp, ct, kexc, hp& +& , ht, q, l) IMPLICIT NONE ! fixed NN output size - REAL(sp), DIMENSION(5), INTENT(IN) :: f_q + REAL(sp), DIMENSION(5), INTENT(IN) :: fq REAL(sp), INTENT(IN) :: pn, en, cp, ct, kexc REAL(sp), INTENT(INOUT) :: hp, ht, q REAL(sp), INTENT(OUT) :: l @@ -14697,20 +14187,20 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE(f_q, pn, en, cp, ct, kexc, & ! dt = 1._sp/real(n_subtimesteps, sp) dt = 1._sp !do i = 1, n_subtimesteps - fhp = (1._sp+f_q(1))*pn*(1._sp-hp**2) - (1._sp+f_q(2))*en*hp*(2._sp-& -& hp) + fhp = (1._sp+fq(1))*pn*(1._sp-hp**2) - (1._sp+fq(2))*en*hp*(2._sp-hp& +& ) hp = hp + dt*fhp*inv_cp IF (hp .LE. 0._sp) hp = 1.e-6_sp IF (hp .GE. 1._sp) hp = 1._sp - 1.e-6_sp - fht = 0.9_sp*(1._sp-f_q(3)**2)*(1._sp+f_q(1))*pn*hp**2 + (1._sp+f_q(& -& 4))*kexc*ht**3.5_sp - (1._sp+f_q(5))*ct*ht**5 + fht = 0.9_sp*(1._sp-fq(3)**2)*(1._sp+fq(1))*pn*hp**2 + (1._sp+fq(4))& +& *kexc*ht**3.5_sp - (1._sp+fq(5))*ct*ht**5 ht = ht + dt*fht/ct IF (ht .LE. 0._sp) ht = 1.e-6_sp IF (ht .GE. 1._sp) ht = 1._sp - 1.e-6_sp !end do - l = (1._sp+f_q(4))*kexc*ht**3.5_sp - q = (1._sp+f_q(5))*ct*ht**5 + (0.1_sp+0.9_sp*f_q(3)**2)*(1._sp+f_q(1& -& ))*pn*hp**2 + l + l = (1._sp+fq(4))*kexc*ht**3.5_sp + q = (1._sp+fq(5))*ct*ht**5 + (0.1_sp+0.9_sp*fq(3)**2)*(1._sp+fq(1))*& +& pn*hp**2 + l END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE ! Differentiation of gr4_time_step in forward (tangent) mode (with options fixinterface noISIZE context): @@ -14763,11 +14253,11 @@ SUBROUTINE GR4_TIME_STEP_D(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION_D(ac_prcp(k), ac_prcp_d(k), ac_pet(k), & & ac_ci(k), ac_ci_d(k), ac_hi(k), ac_hi_d(k)& & , pn, pn_d, en, en_d) - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k& -& ), beta, ac_hp(k), ac_hp_d(k), pr, pr_d, perc& -& , perc_d) - CALL GR_EXCHANGE_D(ac_kexc(k), ac_kexc_d(k), ac_ht(k), & -& ac_ht_d(k), l, l_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_cp(k), ac_cp_d(k), beta, ac_hp(k& +& ), ac_hp_d(k), pr, pr_d, perc, perc_d) + CALL GR_EXCHANGE_D(0._sp, 0.0_4, ac_kexc(k), ac_kexc_d(k), & +& ac_ht(k), ac_ht_d(k), l, l_d) ELSE pr = 0._sp perc = 0._sp @@ -14833,6 +14323,9 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, returns, & REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd REAL(sp) :: pn_b, en_b, pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 + REAL(sp) :: dummydiff_b1 INTEGER :: branch CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) @@ -14855,9 +14348,9 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) - CALL GR_EXCHANGE(ac_kexc(k), ac_ht(k), l) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) + CALL GR_EXCHANGE(0._sp, ac_kexc(k), ac_ht(k), l) CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) @@ -14908,12 +14401,15 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, returns, & l_b = l_b + prr_b CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN - CALL GR_EXCHANGE_B(ac_kexc(k), ac_kexc_b(k), ac_ht(k), & -& ac_ht_b(k), l, l_b) + CALL GR_EXCHANGE_B(0._sp, dummydiff_b1, ac_kexc(k), & +& ac_kexc_b(k), ac_ht(k), ac_ht_b(k), l, l_b) CALL POPREAL4(ac_hp(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k& -& ), beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc& -& , perc_b) + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k), & +& beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc, & +& perc_b) CALL POPREAL4(ac_hi(k)) CALL POPREAL4(pn) CALL POPREAL4(en) @@ -14961,9 +14457,9 @@ SUBROUTINE GR4_TIME_STEP(setup, mesh, input_data, options, returns, & IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) - CALL GR_EXCHANGE(ac_kexc(k), ac_ht(k), l) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) + CALL GR_EXCHANGE(0._sp, ac_kexc(k), ac_ht(k), l) ELSE pr = 0._sp perc = 0._sp @@ -15037,6 +14533,8 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_D(setup, mesh, input_data, options, & REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_d, pn_d, en_d INTEGER :: row, col, k, time_step_returns REAL(sp) :: beta, pr, perc, l, prr, prd, qr, qd + REAL(sp) :: pr_d, perc_d, l_d, prr_d, prd_d, qr_d, qd_d + INTRINSIC MAX REAL(sp) :: temp CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) @@ -15068,25 +14566,24 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_D(setup, mesh, input_data, options, & END DO END DO output_layer_d = 0.0_4 - input_layer_d = 0.0_4 ! Forward MLP without OPENMP DO col=1,mesh%ncol DO row=1,mesh%nrow IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - input_layer_d(1) = ac_hp_d(k) - input_layer(1) = ac_hp(k) - input_layer_d(2) = ac_ht_d(k) - input_layer(2) = ac_ht(k) - input_layer_d(3) = pn_d(k) - input_layer(3) = pn(k) - input_layer_d(4) = en_d(k) - input_layer(4) = en(k) - CALL FORWARD_MLP_D(weight_1, weight_1_d, bias_1, bias_1_d, & -& weight_2, weight_2_d, bias_2, bias_2_d, & -& input_layer, input_layer_d, output_layer(:, k), & -& output_layer_d(:, k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + input_layer_d(:) = (/ac_hp_d(k), ac_ht_d(k), pn_d(k), en_d(k& +& )/) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP_D(weight_1, weight_1_d, bias_1, bias_1_d, & +& weight_2, weight_2_d, bias_2, bias_2_d, & +& input_layer, input_layer_d, output_layer(:, k)& +& , output_layer_d(:, k)) + ELSE + output_layer_d(:, k) = 0.0_4 + output_layer(:, k) = 0._sp + END IF END IF END DO END DO @@ -15096,16 +14593,42 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_D(setup, mesh, input_data, options, & IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - CALL GR_PRODUCTION_TRANSFER_MLP_ALG_D(output_layer(:, k), & -& output_layer_d(:, k), pn(k), & -& pn_d(k), en(k), en_d(k), ac_cp& -& (k), ac_cp_d(k), beta, ac_ct(k& -& ), ac_ct_d(k), ac_kexc(k), & -& ac_kexc_d(k), 5._sp, ac_prcp(k& -& ), ac_hp(k), ac_hp_d(k), ac_ht& -& (k), ac_ht_d(k), ac_qt(k), & -& ac_qt_d(k), pr, perc, l, prr, & -& prd, qr, qd) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL GR_PRODUCTION_D(output_layer(1, k), output_layer_d(1, k& +& ), output_layer(2, k), output_layer_d(2, k), & +& pn(k), pn_d(k), en(k), en_d(k), ac_cp(k), & +& ac_cp_d(k), beta, ac_hp(k), ac_hp_d(k), pr, & +& pr_d, perc, perc_d) + CALL GR_EXCHANGE_D(output_layer(4, k), output_layer_d(4, k)& +& , ac_kexc(k), ac_kexc_d(k), ac_ht(k), ac_ht_d(k& +& ), l, l_d) + ELSE + pr = 0._sp + perc = 0._sp + l = 0._sp + l_d = 0.0_4 + perc_d = 0.0_4 + pr_d = 0.0_4 + END IF + temp = -(output_layer(3, k)*output_layer(3, k)) + 1._sp + prr_d = 0.9_sp*(temp*(pr_d+perc_d)-(pr+perc)*2*output_layer(3& +& , k)*output_layer_d(3, k)) + l_d + prr = 0.9_sp*(temp*(pr+perc)) + l + temp = 0.9_sp*(output_layer(3, k)*output_layer(3, k)) + 0.1_sp + prd_d = (pr+perc)*0.9_sp*2*output_layer(3, k)*output_layer_d(3& +& , k) + temp*(pr_d+perc_d) + prd = temp*(pr+perc) + CALL GR_TRANSFER_D(5._sp, ac_prcp(k), prr, prr_d, ac_ct(k), & +& ac_ct_d(k), ac_ht(k), ac_ht_d(k), qr, qr_d) + IF (0._sp .LT. prd + l) THEN + qd_d = prd_d + l_d + qd = prd + l + ELSE + qd = 0._sp + qd_d = 0.0_4 + END IF + ac_qt_d(k) = qr_d + qd_d + ac_qt(k) = qr + qd ! Transform from mm/dt to m3/s temp = 1e-3_sp*mesh%dx(row, col)*mesh%dy(row, col) ac_qt_d(k) = temp*ac_qt_d(k)/setup%dt @@ -15167,6 +14690,9 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_B(setup, mesh, input_data, options, & REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_b, pn_b, en_b INTEGER :: row, col, k, time_step_returns REAL(sp) :: beta, pr, perc, l, prr, prd, qr, qd + REAL(sp) :: pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b + INTRINSIC MAX + REAL(sp) :: temp_b INTEGER :: branch CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) @@ -15201,20 +14727,19 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_B(setup, mesh, input_data, options, & DO row=1,mesh%nrow IF (mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0) THEN - CALL PUSHCONTROL1B(0) + CALL PUSHCONTROL2B(0) ELSE k = mesh%rowcol_to_ind_ac(row, col) - CALL PUSHREAL4(input_layer(1)) - input_layer(1) = ac_hp(k) - CALL PUSHREAL4(input_layer(2)) - input_layer(2) = ac_ht(k) - CALL PUSHREAL4(input_layer(3)) - input_layer(3) = pn(k) - CALL PUSHREAL4(input_layer(4)) - input_layer(4) = en(k) - CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & -& input_layer, output_layer(:, k)) - CALL PUSHCONTROL1B(1) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL PUSHREAL4ARRAY(input_layer, 4) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & +& input_layer, output_layer(:, k)) + CALL PUSHCONTROL2B(2) + ELSE + output_layer(:, k) = 0._sp + CALL PUSHCONTROL2B(1) + END IF END IF END DO END DO @@ -15225,16 +14750,37 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_B(setup, mesh, input_data, options, & & local_active_cell(row, col) .EQ. 0) THEN CALL PUSHCONTROL1B(0) ELSE + CALL PUSHINTEGER4(k) k = mesh%rowcol_to_ind_ac(row, col) - CALL PUSHREAL4(ac_qt(k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL PUSHREAL4(perc) + CALL PUSHREAL4(pr) + CALL PUSHREAL4(ac_hp(k)) + CALL GR_PRODUCTION(output_layer(1, k), output_layer(2, k), & +& pn(k), en(k), ac_cp(k), beta, ac_hp(k), pr, & +& perc) + CALL GR_EXCHANGE(output_layer(4, k), ac_kexc(k), ac_ht(k), l& +& ) + CALL PUSHCONTROL1B(1) + ELSE + CALL PUSHREAL4(pr) + pr = 0._sp + CALL PUSHREAL4(perc) + perc = 0._sp + l = 0._sp + CALL PUSHCONTROL1B(0) + END IF + CALL PUSHREAL4(prr) + prr = 0.9_sp*(1._sp-output_layer(3, k)**2)*(pr+perc) + l + prd = (0.1_sp+0.9_sp*output_layer(3, k)**2)*(pr+perc) CALL PUSHREAL4(ac_ht(k)) - CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION_TRANSFER_MLP_ALG(output_layer(:, k), pn(k)& -& , en(k), ac_cp(k), beta, ac_ct(k& -& ), ac_kexc(k), 5._sp, ac_prcp(k)& -& , ac_hp(k), ac_ht(k), ac_qt(k), & -& pr, perc, l, prr, prd, qr, qd) -! Transform from mm/dt to m3/s + CALL GR_TRANSFER(5._sp, ac_prcp(k), prr, ac_ct(k), ac_ht(k), & +& qr) + IF (0._sp .LT. prd + l) THEN + CALL PUSHCONTROL1B(0) + ELSE + CALL PUSHCONTROL1B(1) + END IF CALL PUSHCONTROL1B(1) END IF END DO @@ -15246,49 +14792,74 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_B(setup, mesh, input_data, options, & DO row=mesh%nrow,1,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN - k = mesh%rowcol_to_ind_ac(row, col) ac_qt_b(k) = mesh%dx(row, col)*1e-3_sp*mesh%dy(row, col)*& & ac_qt_b(k)/setup%dt - CALL POPREAL4(ac_hp(k)) - CALL POPREAL4(ac_ht(k)) - CALL POPREAL4(ac_qt(k)) - CALL GR_PRODUCTION_TRANSFER_MLP_ALG_B(output_layer(:, k), & -& output_layer_b(:, k), pn(k), & -& pn_b(k), en(k), en_b(k), ac_cp& -& (k), ac_cp_b(k), beta, ac_ct(k& -& ), ac_ct_b(k), ac_kexc(k), & -& ac_kexc_b(k), 5._sp, ac_prcp(k& -& ), ac_hp(k), ac_hp_b(k), ac_ht& -& (k), ac_ht_b(k), ac_qt(k), & -& ac_qt_b(k), pr, perc, l, prr, & -& prd, qr, qd) + qr_b = ac_qt_b(k) + qd_b = ac_qt_b(k) ac_qt_b(k) = 0.0_4 + CALL POPCONTROL1B(branch) + IF (branch .EQ. 0) THEN + prd_b = qd_b + l_b = qd_b + ELSE + l_b = 0.0_4 + prd_b = 0.0_4 + END IF + CALL POPREAL4(ac_ht(k)) + CALL GR_TRANSFER_B(5._sp, ac_prcp(k), prr, prr_b, ac_ct(k), & +& ac_ct_b(k), ac_ht(k), ac_ht_b(k), qr, qr_b) + output_layer_b(3, k) = output_layer_b(3, k) + 2*output_layer(3& +& , k)*0.9_sp*(pr+perc)*prd_b - 2*output_layer(3, k)*(pr+perc)& +& *0.9_sp*prr_b + temp_b = (0.9_sp*output_layer(3, k)**2+0.1_sp)*prd_b + pr_b = temp_b + perc_b = temp_b + CALL POPREAL4(prr) + temp_b = (1._sp-output_layer(3, k)**2)*0.9_sp*prr_b + l_b = l_b + prr_b + pr_b = pr_b + temp_b + perc_b = perc_b + temp_b + CALL POPCONTROL1B(branch) + IF (branch .EQ. 0) THEN + CALL POPREAL4(perc) + CALL POPREAL4(pr) + ELSE + CALL GR_EXCHANGE_B(output_layer(4, k), output_layer_b(4, k)& +& , ac_kexc(k), ac_kexc_b(k), ac_ht(k), ac_ht_b(k& +& ), l, l_b) + CALL POPREAL4(ac_hp(k)) + CALL POPREAL4(pr) + CALL POPREAL4(perc) + CALL GR_PRODUCTION_B(output_layer(1, k), output_layer_b(1, k& +& ), output_layer(2, k), output_layer_b(2, k), & +& pn(k), pn_b(k), en(k), en_b(k), ac_cp(k), & +& ac_cp_b(k), beta, ac_hp(k), ac_hp_b(k), pr, & +& pr_b, perc, perc_b) + END IF + CALL POPINTEGER4(k) END IF END DO END DO - input_layer_b = 0.0_4 DO col=mesh%ncol,1,-1 DO row=mesh%nrow,1,-1 - CALL POPCONTROL1B(branch) + CALL POPCONTROL2B(branch) IF (branch .NE. 0) THEN - k = mesh%rowcol_to_ind_ac(row, col) - CALL FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & -& weight_2, weight_2_b, bias_2, bias_2_b, & -& input_layer, input_layer_b, output_layer(:, k), & -& output_layer_b(:, k)) - output_layer_b(:, k) = 0.0_4 - CALL POPREAL4(input_layer(4)) - en_b(k) = en_b(k) + input_layer_b(4) - input_layer_b(4) = 0.0_4 - CALL POPREAL4(input_layer(3)) - pn_b(k) = pn_b(k) + input_layer_b(3) - input_layer_b(3) = 0.0_4 - CALL POPREAL4(input_layer(2)) - ac_ht_b(k) = ac_ht_b(k) + input_layer_b(2) - input_layer_b(2) = 0.0_4 - CALL POPREAL4(input_layer(1)) - ac_hp_b(k) = ac_hp_b(k) + input_layer_b(1) - input_layer_b(1) = 0.0_4 + IF (branch .EQ. 1) THEN + k = mesh%rowcol_to_ind_ac(row, col) + output_layer_b(:, k) = 0.0_4 + ELSE + k = mesh%rowcol_to_ind_ac(row, col) + CALL FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & +& weight_2, weight_2_b, bias_2, bias_2_b, & +& input_layer, input_layer_b, output_layer(:, k)& +& , output_layer_b(:, k)) + output_layer_b(:, k) = 0.0_4 + CALL POPREAL4ARRAY(input_layer, 4) + ac_hp_b(k) = ac_hp_b(k) + input_layer_b(1) + ac_ht_b(k) = ac_ht_b(k) + input_layer_b(2) + pn_b(k) = pn_b(k) + input_layer_b(3) + en_b(k) = en_b(k) + input_layer_b(4) + END IF END IF END DO END DO @@ -15344,6 +14915,7 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP(setup, mesh, input_data, options, & REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet, pn, en INTEGER :: row, col, k, time_step_returns REAL(sp) :: beta, pr, perc, l, prr, prd, qr, qd + INTRINSIC MAX CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& @@ -15373,12 +14945,13 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP(setup, mesh, input_data, options, & IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - input_layer(1) = ac_hp(k) - input_layer(2) = ac_ht(k) - input_layer(3) = pn(k) - input_layer(4) = en(k) - CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & -& input_layer, output_layer(:, k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & +& input_layer, output_layer(:, k)) + ELSE + output_layer(:, k) = 0._sp + END IF END IF END DO END DO @@ -15388,11 +14961,27 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP(setup, mesh, input_data, options, & IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - CALL GR_PRODUCTION_TRANSFER_MLP_ALG(output_layer(:, k), pn(k)& -& , en(k), ac_cp(k), beta, ac_ct(k& -& ), ac_kexc(k), 5._sp, ac_prcp(k)& -& , ac_hp(k), ac_ht(k), ac_qt(k), & -& pr, perc, l, prr, prd, qr, qd) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL GR_PRODUCTION(output_layer(1, k), output_layer(2, k), & +& pn(k), en(k), ac_cp(k), beta, ac_hp(k), pr, & +& perc) + CALL GR_EXCHANGE(output_layer(4, k), ac_kexc(k), ac_ht(k), l& +& ) + ELSE + pr = 0._sp + perc = 0._sp + l = 0._sp + END IF + prr = 0.9_sp*(1._sp-output_layer(3, k)**2)*(pr+perc) + l + prd = (0.1_sp+0.9_sp*output_layer(3, k)**2)*(pr+perc) + CALL GR_TRANSFER(5._sp, ac_prcp(k), prr, ac_ct(k), ac_ht(k), & +& qr) + IF (0._sp .LT. prd + l) THEN + qd = prd + l + ELSE + qd = 0._sp + END IF + ac_qt(k) = qr + qd ! Transform from mm/dt to m3/s ac_qt(k) = ac_qt(k)*1e-3_sp*mesh%dx(row, col)*mesh%dy(row, col& & )/setup%dt @@ -15738,25 +15327,24 @@ SUBROUTINE GR4_MLP_ODE_TIME_STEP_D(setup, mesh, input_data, options, & END DO END DO output_layer_d = 0.0_4 - input_layer_d = 0.0_4 ! Forward MLP without OPENMP DO col=1,mesh%ncol DO row=1,mesh%nrow IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - input_layer_d(1) = ac_hp_d(k) - input_layer(1) = ac_hp(k) - input_layer_d(2) = ac_ht_d(k) - input_layer(2) = ac_ht(k) - input_layer_d(3) = pn_d(k) - input_layer(3) = pn(k) - input_layer_d(4) = en_d(k) - input_layer(4) = en(k) - CALL FORWARD_MLP_D(weight_1, weight_1_d, bias_1, bias_1_d, & -& weight_2, weight_2_d, bias_2, bias_2_d, & -& input_layer, input_layer_d, output_layer(:, k), & -& output_layer_d(:, k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + input_layer_d(:) = (/ac_hp_d(k), ac_ht_d(k), pn_d(k), en_d(k& +& )/) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP_D(weight_1, weight_1_d, bias_1, bias_1_d, & +& weight_2, weight_2_d, bias_2, bias_2_d, & +& input_layer, input_layer_d, output_layer(:, k)& +& , output_layer_d(:, k)) + ELSE + output_layer_d(:, k) = 0.0_4 + output_layer(:, k) = 0._sp + END IF END IF END DO END DO @@ -15867,20 +15455,19 @@ SUBROUTINE GR4_MLP_ODE_TIME_STEP_B(setup, mesh, input_data, options, & DO row=1,mesh%nrow IF (mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0) THEN - CALL PUSHCONTROL1B(0) + CALL PUSHCONTROL2B(0) ELSE k = mesh%rowcol_to_ind_ac(row, col) - CALL PUSHREAL4(input_layer(1)) - input_layer(1) = ac_hp(k) - CALL PUSHREAL4(input_layer(2)) - input_layer(2) = ac_ht(k) - CALL PUSHREAL4(input_layer(3)) - input_layer(3) = pn(k) - CALL PUSHREAL4(input_layer(4)) - input_layer(4) = en(k) - CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & -& input_layer, output_layer(:, k)) - CALL PUSHCONTROL1B(1) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL PUSHREAL4ARRAY(input_layer, 4) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & +& input_layer, output_layer(:, k)) + CALL PUSHCONTROL2B(2) + ELSE + output_layer(:, k) = 0._sp + CALL PUSHCONTROL2B(1) + END IF END IF END DO END DO @@ -15929,29 +15516,26 @@ SUBROUTINE GR4_MLP_ODE_TIME_STEP_B(setup, mesh, input_data, options, & END IF END DO END DO - input_layer_b = 0.0_4 DO col=mesh%ncol,1,-1 DO row=mesh%nrow,1,-1 - CALL POPCONTROL1B(branch) + CALL POPCONTROL2B(branch) IF (branch .NE. 0) THEN - k = mesh%rowcol_to_ind_ac(row, col) - CALL FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & -& weight_2, weight_2_b, bias_2, bias_2_b, & -& input_layer, input_layer_b, output_layer(:, k), & -& output_layer_b(:, k)) - output_layer_b(:, k) = 0.0_4 - CALL POPREAL4(input_layer(4)) - en_b(k) = en_b(k) + input_layer_b(4) - input_layer_b(4) = 0.0_4 - CALL POPREAL4(input_layer(3)) - pn_b(k) = pn_b(k) + input_layer_b(3) - input_layer_b(3) = 0.0_4 - CALL POPREAL4(input_layer(2)) - ac_ht_b(k) = ac_ht_b(k) + input_layer_b(2) - input_layer_b(2) = 0.0_4 - CALL POPREAL4(input_layer(1)) - ac_hp_b(k) = ac_hp_b(k) + input_layer_b(1) - input_layer_b(1) = 0.0_4 + IF (branch .EQ. 1) THEN + k = mesh%rowcol_to_ind_ac(row, col) + output_layer_b(:, k) = 0.0_4 + ELSE + k = mesh%rowcol_to_ind_ac(row, col) + CALL FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & +& weight_2, weight_2_b, bias_2, bias_2_b, & +& input_layer, input_layer_b, output_layer(:, k)& +& , output_layer_b(:, k)) + output_layer_b(:, k) = 0.0_4 + CALL POPREAL4ARRAY(input_layer, 4) + ac_hp_b(k) = ac_hp_b(k) + input_layer_b(1) + ac_ht_b(k) = ac_ht_b(k) + input_layer_b(2) + pn_b(k) = pn_b(k) + input_layer_b(3) + en_b(k) = en_b(k) + input_layer_b(4) + END IF END IF END DO END DO @@ -16034,12 +15618,13 @@ SUBROUTINE GR4_MLP_ODE_TIME_STEP(setup, mesh, input_data, options, & IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - input_layer(1) = ac_hp(k) - input_layer(2) = ac_ht(k) - input_layer(3) = pn(k) - input_layer(4) = en(k) - CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & -& input_layer, output_layer(:, k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & +& input_layer, output_layer(:, k)) + ELSE + output_layer(:, k) = 0._sp + END IF END IF END DO END DO @@ -16111,9 +15696,9 @@ SUBROUTINE GR5_TIME_STEP_D(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION_D(ac_prcp(k), ac_prcp_d(k), ac_pet(k), & & ac_ci(k), ac_ci_d(k), ac_hi(k), ac_hi_d(k)& & , pn, pn_d, en, en_d) - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k& -& ), beta, ac_hp(k), ac_hp_d(k), pr, pr_d, perc& -& , perc_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_cp(k), ac_cp_d(k), beta, ac_hp(k& +& ), ac_hp_d(k), pr, pr_d, perc, perc_d) CALL GR_THRESHOLD_EXCHANGE_D(ac_kexc(k), ac_kexc_d(k), & & ac_aexc(k), ac_aexc_d(k), ac_ht(k), & & ac_ht_d(k), l, l_d) @@ -16182,6 +15767,8 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, returns, & REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd REAL(sp) :: pn_b, en_b, pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 INTEGER :: branch CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) @@ -16204,8 +15791,8 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)& & , l) CALL PUSHCONTROL1B(1) @@ -16262,9 +15849,12 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, returns, & & ac_aexc(k), ac_aexc_b(k), ac_ht(k), & & ac_ht_b(k), l, l_b) CALL POPREAL4(ac_hp(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k& -& ), beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc& -& , perc_b) + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k), & +& beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc, & +& perc_b) CALL POPREAL4(ac_hi(k)) CALL POPREAL4(pn) CALL POPREAL4(en) @@ -16312,8 +15902,8 @@ SUBROUTINE GR5_TIME_STEP(setup, mesh, input_data, options, returns, & IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)& & , l) ELSE @@ -16392,9 +15982,9 @@ SUBROUTINE GR6_TIME_STEP_D(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION_D(ac_prcp(k), ac_prcp_d(k), ac_pet(k), & & ac_ci(k), ac_ci_d(k), ac_hi(k), ac_hi_d(k)& & , pn, pn_d, en, en_d) - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k& -& ), beta, ac_hp(k), ac_hp_d(k), pr, pr_d, perc& -& , perc_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_cp(k), ac_cp_d(k), beta, ac_hp(k& +& ), ac_hp_d(k), pr, pr_d, perc, perc_d) CALL GR_THRESHOLD_EXCHANGE_D(ac_kexc(k), ac_kexc_d(k), & & ac_aexc(k), ac_aexc_d(k), ac_ht(k), & & ac_ht_d(k), l, l_d) @@ -16470,6 +16060,8 @@ SUBROUTINE GR6_TIME_STEP_B(setup, mesh, input_data, options, returns, & REAL(sp) :: pn_b, en_b, pr_b, perc_b, l_b, prr_b, pre_b, prd_b, qr_b& & , qd_b, qe_b INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 REAL(sp) :: temp_b INTEGER :: branch CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& @@ -16493,8 +16085,8 @@ SUBROUTINE GR6_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)& & , l) CALL PUSHCONTROL1B(1) @@ -16564,9 +16156,12 @@ SUBROUTINE GR6_TIME_STEP_B(setup, mesh, input_data, options, returns, & & ac_aexc(k), ac_aexc_b(k), ac_ht(k), & & ac_ht_b(k), l, l_b) CALL POPREAL4(ac_hp(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k& -& ), beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc& -& , perc_b) + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k), & +& beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc, & +& perc_b) CALL POPREAL4(ac_hi(k)) CALL POPREAL4(pn) CALL POPREAL4(en) @@ -16615,8 +16210,8 @@ SUBROUTINE GR6_TIME_STEP(setup, mesh, input_data, options, returns, & IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)& & , l) ELSE @@ -16702,9 +16297,9 @@ SUBROUTINE GRD_TIME_STEP_D(setup, mesh, input_data, options, returns, & END IF en_d = -ei_d en = ac_pet(k) - ei - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k& -& ), 1000._sp, ac_hp(k), ac_hp_d(k), pr, pr_d, & -& perc, perc_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_cp(k), ac_cp_d(k), 1000._sp, & +& ac_hp(k), ac_hp_d(k), pr, pr_d, perc, perc_d) ELSE pr = 0._sp perc = 0._sp @@ -16756,6 +16351,8 @@ SUBROUTINE GRD_TIME_STEP_B(setup, mesh, input_data, options, returns, & REAL(sp) :: ei_b, pn_b, en_b, pr_b, perc_b, prr_b, qr_b INTRINSIC MIN INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 INTEGER :: branch CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) @@ -16789,8 +16386,8 @@ SUBROUTINE GRD_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL PUSHREAL4(en) en = ac_pet(k) - ei CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION(pn, en, ac_cp(k), 1000._sp, ac_hp(k), pr& -& , perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), 1000._sp& +& , ac_hp(k), pr, perc) CALL PUSHCONTROL1B(0) ELSE CALL PUSHCONTROL1B(1) @@ -16826,8 +16423,11 @@ SUBROUTINE GRD_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN CALL POPREAL4(ac_hp(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k& -& ), 1000._sp, ac_hp(k), ac_hp_b(k), pr, pr_b, & + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k), & +& 1000._sp, ac_hp(k), ac_hp_b(k), pr, pr_b, & & perc, perc_b) CALL POPREAL4(en) ei_b = -en_b @@ -16888,8 +16488,8 @@ SUBROUTINE GRD_TIME_STEP(setup, mesh, input_data, options, returns, & pn = 0._sp END IF en = ac_pet(k) - ei - CALL GR_PRODUCTION(pn, en, ac_cp(k), 1000._sp, ac_hp(k), pr& -& , perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), 1000._sp& +& , ac_hp(k), pr, perc) ELSE pr = 0._sp perc = 0._sp @@ -16967,9 +16567,9 @@ SUBROUTINE LOIEAU_TIME_STEP_D(setup, mesh, input_data, options, & END IF en_d = -ei_d en = ac_pet(k) - ei - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_ca(k), ac_ca_d(k& -& ), beta, ac_ha(k), ac_ha_d(k), pr, pr_d, perc& -& , perc_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_ca(k), ac_ca_d(k), beta, ac_ha(k& +& ), ac_ha_d(k), pr, pr_d, perc, perc_d) ELSE pr = 0._sp perc = 0._sp @@ -17030,6 +16630,8 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, & REAL(sp) :: ei_b, pn_b, en_b, pr_b, perc_b, prr_b, prd_b, qr_b, qd_b INTRINSIC MIN INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 INTEGER :: branch CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) @@ -17065,8 +16667,8 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, & CALL PUSHREAL4(en) en = ac_pet(k) - ei CALL PUSHREAL4(ac_ha(k)) - CALL GR_PRODUCTION(pn, en, ac_ca(k), beta, ac_ha(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_ca(k), beta, & +& ac_ha(k), pr, perc) CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) @@ -17123,9 +16725,12 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, & CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN CALL POPREAL4(ac_ha(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_ca(k), ac_ca_b(k& -& ), beta, ac_ha(k), ac_ha_b(k), pr, pr_b, perc& -& , perc_b) + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_ca(k), ac_ca_b(k), & +& beta, ac_ha(k), ac_ha_b(k), pr, pr_b, perc, & +& perc_b) CALL POPREAL4(en) ei_b = -en_b CALL POPCONTROL1B(branch) @@ -17187,8 +16792,8 @@ SUBROUTINE LOIEAU_TIME_STEP(setup, mesh, input_data, options, returns& pn = 0._sp END IF en = ac_pet(k) - ei - CALL GR_PRODUCTION(pn, en, ac_ca(k), beta, ac_ha(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_ca(k), beta, & +& ac_ha(k), pr, perc) ELSE pr = 0._sp perc = 0._sp diff --git a/smash/fcore/forward/forward_openmp_db.f90 b/smash/fcore/forward/forward_openmp_db.f90 index 0cca0e94..1a5fd51a 100644 --- a/smash/fcore/forward/forward_openmp_db.f90 +++ b/smash/fcore/forward/forward_openmp_db.f90 @@ -12587,7 +12587,7 @@ SUBROUTINE DOT_PRODUCT_2D_1D_D(a, a_d, x, x_d, b, b_d) END SUBROUTINE DOT_PRODUCT_2D_1D_D ! Differentiation of dot_product_2d_1d in reverse (adjoint) mode (with options fixinterface noISIZE context OpenMP): -! gradient of useful results: x a b +! gradient of useful results: a b ! with respect to varying inputs: x a SUBROUTINE DOT_PRODUCT_2D_1D_B(a, a_b, x, x_b, b, b_b) IMPLICIT NONE @@ -12608,6 +12608,7 @@ SUBROUTINE DOT_PRODUCT_2D_1D_B(a, a_b, x, x_b, b, b_b) CALL PUSHINTEGER4(i - 1) END DO ad_to0 = j - 1 + x_b = 0.0_4 DO j=ad_to0,1,-1 CALL POPINTEGER4(ad_to) DO i=ad_to,1,-1 @@ -12698,7 +12699,7 @@ END SUBROUTINE FORWARD_MLP_D ! Differentiation of forward_mlp in reverse (adjoint) mode (with options fixinterface noISIZE context OpenMP): ! gradient of useful results: output_layer bias_1 bias_2 -! input_layer weight_1 weight_2 +! weight_1 weight_2 ! with respect to varying inputs: bias_1 bias_2 input_layer weight_1 ! weight_2 SUBROUTINE FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & @@ -12750,7 +12751,6 @@ SUBROUTINE FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & output_layer_b(i) = temp_b bias_2_b(i) = bias_2_b(i) + temp_b END DO - inter_layer_b = 0.0_4 CALL DOT_PRODUCT_2D_1D_B(weight_2, weight_2_b, inter_layer, & & inter_layer_b, output_layer, output_layer_b) CALL POPINTEGER4(ad_to) @@ -12805,7 +12805,6 @@ END MODULE MD_NEURAL_NETWORK_DIFF !% - gr_exchange !% - gr_threshold_exchange !% - gr_transfer -!% - gr_production_transfer_mlp_alg !% - gr_production_transfer_ode !% - gr_production_transfer_mlp_ode !% - gr4_time_step @@ -12966,12 +12965,12 @@ END SUBROUTINE GR_INTERCEPTION ! Differentiation of gr_production in forward (tangent) mode (with options fixinterface noISIZE context OpenMP): ! variations of useful results: hp perc pr -! with respect to varying inputs: hp en cp pn - SUBROUTINE GR_PRODUCTION_D(pn, pn_d, en, en_d, cp, cp_d, beta, hp, & -& hp_d, pr, pr_d, perc, perc_d) +! with respect to varying inputs: fq_ps hp en fq_es cp pn + SUBROUTINE GR_PRODUCTION_D(fq_ps, fq_ps_d, fq_es, fq_es_d, pn, pn_d, & +& en, en_d, cp, cp_d, beta, hp, hp_d, pr, pr_d, perc, perc_d) IMPLICIT NONE - REAL(sp), INTENT(IN) :: pn, en, cp, beta - REAL(sp), INTENT(IN) :: pn_d, en_d, cp_d + REAL(sp), INTENT(IN) :: fq_ps, fq_es, pn, en, cp, beta + REAL(sp), INTENT(IN) :: fq_ps_d, fq_es_d, pn_d, en_d, cp_d REAL(sp), INTENT(INOUT) :: hp REAL(sp), INTENT(INOUT) :: hp_d REAL(sp), INTENT(OUT) :: pr, perc @@ -12998,6 +12997,8 @@ SUBROUTINE GR_PRODUCTION_D(pn, pn_d, en, en_d, cp, cp_d, beta, hp, & & inv_cp)**2)*(inv_cp*pn_d+pn*inv_cp_d)-temp2*(temp*hp_d+hp*(1.0-& & TANH(pn*inv_cp)**2)*(inv_cp*pn_d+pn*inv_cp_d)))/(hp*temp+1._sp) ps = temp2 + ps_d = ps*fq_ps_d + (fq_ps+1._sp)*ps_d + ps = (1._sp+fq_ps)*ps temp2 = TANH(en*inv_cp) temp1 = TANH(en*inv_cp) temp0 = hp*cp*(-hp+2._sp) @@ -13007,6 +13008,8 @@ SUBROUTINE GR_PRODUCTION_D(pn, pn_d, en, en_d, cp, cp_d, beta, hp, & & 1.0-TANH(en*inv_cp)**2)*(inv_cp*en_d+en*inv_cp_d)-temp2*hp_d))/((& & 1._sp-hp)*temp2+1._sp) es = temp + es_d = es*fq_es_d + (fq_es+1._sp)*es_d + es = (1._sp+fq_es)*es hp_imd_d = hp_d + inv_cp*(ps_d-es_d) + (ps-es)*inv_cp_d hp_imd = hp + (ps-es)*inv_cp IF (pn .GT. 0) THEN @@ -13026,13 +13029,14 @@ SUBROUTINE GR_PRODUCTION_D(pn, pn_d, en, en_d, cp, cp_d, beta, hp, & END SUBROUTINE GR_PRODUCTION_D ! Differentiation of gr_production in reverse (adjoint) mode (with options fixinterface noISIZE context OpenMP): -! gradient of useful results: hp cp perc pr -! with respect to varying inputs: hp en cp pn - SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & -& hp_b, pr, pr_b, perc, perc_b) - IMPLICIT NONE - REAL(sp), INTENT(IN) :: pn, en, cp, beta - REAL(sp) :: pn_b, en_b, cp_b +! gradient of useful results: fq_ps hp en fq_es cp pn perc +! pr +! with respect to varying inputs: fq_ps hp en fq_es cp pn + SUBROUTINE GR_PRODUCTION_B(fq_ps, fq_ps_b, fq_es, fq_es_b, pn, pn_b, & +& en, en_b, cp, cp_b, beta, hp, hp_b, pr, pr_b, perc, perc_b) + IMPLICIT NONE + REAL(sp), INTENT(IN) :: fq_ps, fq_es, pn, en, cp, beta + REAL(sp) :: fq_ps_b, fq_es_b, pn_b, en_b, cp_b REAL(sp), INTENT(INOUT) :: hp REAL(sp), INTENT(INOUT) :: hp_b REAL(sp) :: pr, perc @@ -13060,8 +13064,12 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & INTEGER :: branch inv_cp = 1._sp/cp ps = cp*(1._sp-hp*hp)*TANH(pn*inv_cp)/(1._sp+hp*TANH(pn*inv_cp)) + CALL PUSHREAL4(ps) + ps = (1._sp+fq_ps)*ps es = hp*cp*(2._sp-hp)*TANH(en*inv_cp)/(1._sp+(1._sp-hp)*TANH(en*& & inv_cp)) + CALL PUSHREAL4(es) + es = (1._sp+fq_es)*es hp_imd = hp + (ps-es)*inv_cp IF (pn .GT. 0) THEN CALL PUSHCONTROL1B(0) @@ -13070,12 +13078,14 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & END IF pwx1 = 1._sp + (hp_imd/beta)**4 pwr1 = pwx1**(-0.25_sp) + CALL PUSHREAL4(perc) perc = hp_imd*cp*(1._sp-pwr1) pwx1 = 1._sp + (hp_imd/beta)**4 pwr1 = pwx1**(-0.25_sp) inv_cp = 1._sp/cp perc_b = perc_b - inv_cp*hp_b inv_cp_b = -(perc*hp_b) + CALL POPREAL4(perc) !$OMP ATOMIC update cp_b = cp_b + hp_imd*(1._sp-pwr1)*perc_b pwr1_b = -(hp_imd*cp*perc_b) @@ -13084,16 +13094,21 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & & 4 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN - pn_b = pr_b +!$OMP ATOMIC update + pn_b = pn_b + pr_b hp_imd_b = hp_imd_b - cp*pr_b hp_b = cp*pr_b !$OMP ATOMIC update cp_b = cp_b - (hp_imd-hp)*pr_b ELSE hp_b = 0.0_4 - pn_b = 0.0_4 END IF es_b = -(inv_cp*hp_imd_b) + inv_cp_b = inv_cp_b + (ps-es)*hp_imd_b + CALL POPREAL4(es) +!$OMP ATOMIC update + fq_es_b = fq_es_b + es*es_b + es_b = (fq_es+1._sp)*es_b temp4 = TANH(en*inv_cp) temp3 = (-hp+1._sp)*temp4 + 1._sp temp1 = TANH(en*inv_cp) @@ -13107,9 +13122,14 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & ps_b = inv_cp*hp_imd_b temp_b4 = (1.0-TANH(en*inv_cp)**2)*temp0*temp_b3 temp_b5 = (1.0-TANH(en*inv_cp)**2)*(1._sp-hp)*temp_b0 - en_b = inv_cp*temp_b5 + inv_cp*temp_b4 +!$OMP ATOMIC update + en_b = en_b + inv_cp*temp_b5 + inv_cp*temp_b4 !$OMP ATOMIC update cp_b = cp_b + hp*temp_b + CALL POPREAL4(ps) +!$OMP ATOMIC update + fq_ps_b = fq_ps_b + ps*ps_b + ps_b = (fq_ps+1._sp)*ps_b temp = TANH(pn*inv_cp) temp0 = hp*temp + 1._sp temp1 = TANH(pn*inv_cp) @@ -13120,16 +13140,17 @@ SUBROUTINE GR_PRODUCTION_B(pn, pn_b, en, en_b, cp, cp_b, beta, hp, & !$OMP ATOMIC update hp_b = hp_b + temp*temp_b1 - 2*hp*cp*temp1*temp_b temp_b2 = (1.0-TANH(pn*inv_cp)**2)*hp*temp_b1 - inv_cp_b = inv_cp_b + (ps-es)*hp_imd_b + en*temp_b5 + en*temp_b4 + & -& pn*temp_b2 + pn*temp_b0 + inv_cp_b = inv_cp_b + en*temp_b5 + en*temp_b4 + pn*temp_b2 + pn*& +& temp_b0 !$OMP ATOMIC update cp_b = cp_b + (1._sp-hp**2)*temp1*temp_b - inv_cp_b/cp**2 +!$OMP ATOMIC update pn_b = pn_b + inv_cp*temp_b2 + inv_cp*temp_b0 END SUBROUTINE GR_PRODUCTION_B - SUBROUTINE GR_PRODUCTION(pn, en, cp, beta, hp, pr, perc) + SUBROUTINE GR_PRODUCTION(fq_ps, fq_es, pn, en, cp, beta, hp, pr, perc) IMPLICIT NONE - REAL(sp), INTENT(IN) :: pn, en, cp, beta + REAL(sp), INTENT(IN) :: fq_ps, fq_es, pn, en, cp, beta REAL(sp), INTENT(INOUT) :: hp REAL(sp), INTENT(OUT) :: pr, perc REAL(sp) :: inv_cp, ps, es, hp_imd @@ -13139,8 +13160,10 @@ SUBROUTINE GR_PRODUCTION(pn, en, cp, beta, hp, pr, perc) inv_cp = 1._sp/cp pr = 0._sp ps = cp*(1._sp-hp*hp)*TANH(pn*inv_cp)/(1._sp+hp*TANH(pn*inv_cp)) + ps = (1._sp+fq_ps)*ps es = hp*cp*(2._sp-hp)*TANH(en*inv_cp)/(1._sp+(1._sp-hp)*TANH(en*& & inv_cp)) + es = (1._sp+fq_es)*es hp_imd = hp + (ps-es)*inv_cp IF (pn .GT. 0) pr = pn - (hp_imd-hp)*cp pwx1 = 1._sp + (hp_imd/beta)**4 @@ -13151,44 +13174,49 @@ END SUBROUTINE GR_PRODUCTION ! Differentiation of gr_exchange in forward (tangent) mode (with options fixinterface noISIZE context OpenMP): ! variations of useful results: l -! with respect to varying inputs: kexc ht - SUBROUTINE GR_EXCHANGE_D(kexc, kexc_d, ht, ht_d, l, l_d) +! with respect to varying inputs: kexc fq_l ht + SUBROUTINE GR_EXCHANGE_D(fq_l, fq_l_d, kexc, kexc_d, ht, ht_d, l, l_d) IMPLICIT NONE - REAL(sp), INTENT(IN) :: kexc - REAL(sp), INTENT(IN) :: kexc_d + REAL(sp), INTENT(IN) :: fq_l, kexc + REAL(sp), INTENT(IN) :: fq_l_d, kexc_d REAL(sp), INTENT(INOUT) :: ht REAL(sp), INTENT(INOUT) :: ht_d REAL(sp), INTENT(OUT) :: l REAL(sp), INTENT(OUT) :: l_d REAL(sp) :: temp temp = ht**3.5_sp - l_d = temp*kexc_d + kexc*3.5_sp*ht**2.5*ht_d - l = kexc*temp + l_d = temp*(kexc*fq_l_d+(fq_l+1._sp)*kexc_d) + (fq_l+1._sp)*kexc*& +& 3.5_sp*ht**2.5*ht_d + l = (fq_l+1._sp)*kexc*temp END SUBROUTINE GR_EXCHANGE_D ! Differentiation of gr_exchange in reverse (adjoint) mode (with options fixinterface noISIZE context OpenMP): -! gradient of useful results: l kexc ht -! with respect to varying inputs: kexc ht - SUBROUTINE GR_EXCHANGE_B(kexc, kexc_b, ht, ht_b, l, l_b) +! gradient of useful results: l kexc fq_l ht +! with respect to varying inputs: kexc fq_l ht + SUBROUTINE GR_EXCHANGE_B(fq_l, fq_l_b, kexc, kexc_b, ht, ht_b, l, l_b) IMPLICIT NONE - REAL(sp), INTENT(IN) :: kexc - REAL(sp) :: kexc_b + REAL(sp), INTENT(IN) :: fq_l, kexc + REAL(sp) :: fq_l_b, kexc_b REAL(sp), INTENT(INOUT) :: ht REAL(sp), INTENT(INOUT) :: ht_b REAL(sp) :: l REAL(sp) :: l_b + REAL(sp) :: temp_b + temp_b = ht**3.5_sp*l_b !$OMP ATOMIC update - kexc_b = kexc_b + ht**3.5_sp*l_b + ht_b = ht_b + 3.5_sp*ht**2.5*(fq_l+1._sp)*kexc*l_b +!$OMP ATOMIC update + fq_l_b = fq_l_b + kexc*temp_b !$OMP ATOMIC update - ht_b = ht_b + 3.5_sp*ht**2.5*kexc*l_b + kexc_b = kexc_b + (fq_l+1._sp)*temp_b END SUBROUTINE GR_EXCHANGE_B - SUBROUTINE GR_EXCHANGE(kexc, ht, l) + SUBROUTINE GR_EXCHANGE(fq_l, kexc, ht, l) IMPLICIT NONE - REAL(sp), INTENT(IN) :: kexc + REAL(sp), INTENT(IN) :: fq_l, kexc REAL(sp), INTENT(INOUT) :: ht REAL(sp), INTENT(OUT) :: l - l = kexc*ht**3.5_sp + l = (1._sp+fq_l)*kexc*ht**3.5_sp END SUBROUTINE GR_EXCHANGE ! Differentiation of gr_threshold_exchange in forward (tangent) mode (with options fixinterface noISIZE context OpenMP): @@ -13614,607 +13642,50 @@ SUBROUTINE GR_EXPONENTIAL_TRANSFER_B(pre, pre_b, be, be_b, he, he_b, & IF (branch .EQ. 0) THEN ar = he_star/be !$OMP ATOMIC update - be_b = be_b + EXP(ar)*qe_b - ar_b = EXP(ar)*be*qe_b - ELSE IF (branch .EQ. 1) THEN - ar = he_star/be - temp = EXP(ar) - he_star_b = he_star_b + qe_b -!$OMP ATOMIC update - be_b = be_b + qe_b/temp - ar_b = -(EXP(ar)*be*qe_b/temp**2) - ELSE -!$OMP ATOMIC update - be_b = be_b + LOG(arg1)*qe_b - arg1_b = be*qe_b/arg1 - ar = he_star/be - ar_b = EXP(ar)*arg1_b - END IF - he_star_b = he_star_b + ar_b/be -!$OMP ATOMIC update - be_b = be_b - he_star*ar_b/be**2 - he_b = he_star_b - pre_b = he_star_b - END SUBROUTINE GR_EXPONENTIAL_TRANSFER_B - - SUBROUTINE GR_EXPONENTIAL_TRANSFER(pre, be, he, qe) - IMPLICIT NONE - REAL(sp), INTENT(IN) :: pre, be - REAL(sp), INTENT(INOUT) :: he - REAL(sp), INTENT(OUT) :: qe - REAL(sp) :: he_star, ar - INTRINSIC EXP - INTRINSIC LOG - REAL(sp) :: arg1 - he_star = he + pre - ar = he_star/be - IF (ar .LT. -7._sp) THEN - qe = be*EXP(ar) - ELSE IF (ar .GT. 7._sp) THEN - qe = he_star + be/EXP(ar) - ELSE - arg1 = EXP(ar) + 1._sp - qe = be*LOG(arg1) - END IF - he = he_star - qe - END SUBROUTINE GR_EXPONENTIAL_TRANSFER - -! Differentiation of gr_production_transfer_mlp_alg in forward (tangent) mode (with options fixinterface noISIZE context OpenMP) -!: -! variations of useful results: q hp ht -! with respect to varying inputs: kexc hp ht en f_q cp pn ct - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG_D(f_q, f_q_d, pn, pn_d, en, & -& en_d, cp, cp_d, beta, ct, ct_d, kexc, kexc_d, n, prcp, hp, hp_d, ht& -& , ht_d, q, q_d, pr, perc, l, prr, prd, qr, qd) - IMPLICIT NONE -! fixed NN output size - REAL(sp), DIMENSION(4), INTENT(IN) :: f_q - REAL(sp), DIMENSION(4), INTENT(IN) :: f_q_d - REAL(sp), INTENT(IN) :: pn, en, cp, beta, ct, kexc, n, prcp - REAL(sp), INTENT(IN) :: pn_d, en_d, cp_d, ct_d, kexc_d - REAL(sp), INTENT(INOUT) :: hp, ht, q - REAL(sp), INTENT(INOUT) :: hp_d, ht_d, q_d - REAL(sp), INTENT(OUT) :: pr, perc, l, prr, prd, qr, qd - REAL(sp) :: pr_d - REAL(sp) :: inv_cp, ps, es, hp_imd, pr_imd, ht_imd, nm1, d1pnm1 - REAL(sp) :: inv_cp_d, ps_d, es_d, hp_imd_d, pr_imd_d, ht_imd_d - INTRINSIC TANH - INTRINSIC MAX - REAL(sp) :: pwx1 - REAL(sp) :: pwx1_d - REAL(sp) :: pwr1 - REAL(sp) :: pwr1_d - REAL(sp) :: pwy1 - REAL(sp) :: pwy2 - REAL(sp) :: pwr2 - REAL(sp) :: pwr2_d - REAL(sp) :: pwx3 - REAL(sp) :: pwx3_d - REAL(sp) :: pwy3 - REAL(sp) :: pwr3 - REAL(sp) :: pwr3_d - REAL(sp) :: temp - REAL(sp) :: temp0 - REAL(sp) :: temp1 - REAL(sp) :: temp2 - REAL(sp) :: perc_d - REAL(sp) :: qr_d - REAL(sp) :: prd_d - REAL(sp) :: l_d - REAL(sp) :: qd_d - REAL(sp) :: prr_d -!% Production - inv_cp_d = -(cp_d/cp**2) - inv_cp = 1._sp/cp - pr = 0._sp - temp = TANH(pn*inv_cp) - temp0 = TANH(pn*inv_cp) - temp1 = cp*(-(hp*hp)+1._sp) - temp2 = temp1*temp0/(hp*temp+1._sp) - ps_d = (temp0*((1._sp-hp**2)*cp_d-cp*2*hp*hp_d)+temp1*(1.0-TANH(pn*& -& inv_cp)**2)*(inv_cp*pn_d+pn*inv_cp_d)-temp2*(temp*hp_d+hp*(1.0-& -& TANH(pn*inv_cp)**2)*(inv_cp*pn_d+pn*inv_cp_d)))/(hp*temp+1._sp) - ps = temp2 - ps_d = ps*f_q_d(1) + (f_q(1)+1._sp)*ps_d - ps = (1._sp+f_q(1))*ps - temp2 = TANH(en*inv_cp) - temp1 = TANH(en*inv_cp) - temp0 = hp*cp*(-hp+2._sp) - temp = temp0*temp1/((-hp+1._sp)*temp2+1._sp) - es_d = (temp1*((2._sp-hp)*(cp*hp_d+hp*cp_d)-hp*cp*hp_d)+temp0*(1.0-& -& TANH(en*inv_cp)**2)*(inv_cp*en_d+en*inv_cp_d)-temp*((1._sp-hp)*(& -& 1.0-TANH(en*inv_cp)**2)*(inv_cp*en_d+en*inv_cp_d)-temp2*hp_d))/((& -& 1._sp-hp)*temp2+1._sp) - es = temp - es_d = es*f_q_d(2) + (f_q(2)+1._sp)*es_d - es = (1._sp+f_q(2))*es - hp_imd_d = hp_d + inv_cp*(ps_d-es_d) + (ps-es)*inv_cp_d - hp_imd = hp + (ps-es)*inv_cp - IF (pn .GT. 0) THEN - pr_d = pn_d - cp*(hp_imd_d-hp_d) - (hp_imd-hp)*cp_d - pr = pn - (hp_imd-hp)*cp - ELSE - pr_d = 0.0_4 - END IF - pwx1_d = 4*hp_imd**3*hp_imd_d/beta**4 - pwx1 = 1._sp + (hp_imd/beta)**4 - pwr1_d = -(0.25_sp*pwx1**(-1.25)*pwx1_d) - pwr1 = pwx1**(-0.25_sp) - perc_d = (1._sp-pwr1)*(cp*hp_imd_d+hp_imd*cp_d) - hp_imd*cp*pwr1_d - perc = hp_imd*cp*(1._sp-pwr1) - hp_d = hp_imd_d - inv_cp*perc_d - perc*inv_cp_d - hp = hp_imd - perc*inv_cp -!% Transfer - nm1 = n - 1._sp - d1pnm1 = 1._sp/nm1 - temp2 = ht**3.5_sp - l_d = temp2*(kexc*f_q_d(4)+(f_q(4)+1._sp)*kexc_d) + (f_q(4)+1._sp)*& -& kexc*3.5_sp*ht**2.5*ht_d - l = (f_q(4)+1._sp)*kexc*temp2 - prr_d = 0.9_sp*((1._sp-f_q(3)**2)*(pr_d+perc_d)-(pr+perc)*2*f_q(3)*& -& f_q_d(3)) + l_d - prr = 0.9_sp*(1._sp-f_q(3)**2)*(pr+perc) + l - temp2 = 0.9_sp*(f_q(3)*f_q(3)) + 0.1_sp - prd_d = (pr+perc)*0.9_sp*2*f_q(3)*f_q_d(3) + temp2*(pr_d+perc_d) - prd = temp2*(pr+perc) - IF (prcp .LT. 0._sp) THEN - pwx1_d = ct*ht_d + ht*ct_d - pwx1 = ht*ct - pwy1 = -nm1 - IF (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1))) & -& THEN - pwr1_d = 0.0_4 - ELSE - pwr1_d = pwy1*pwx1**(pwy1-1)*pwx1_d - END IF - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - IF (ct .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2))) & -& THEN - pwr2_d = 0.0_4 - ELSE - pwr2_d = pwy2*ct**(pwy2-1)*ct_d - END IF - pwr2 = ct**pwy2 - pwx3_d = pwr1_d - pwr2_d - pwx3 = pwr1 - pwr2 - pwy3 = -d1pnm1 - IF (pwx3 .LE. 0.0 .AND. (pwy3 .EQ. 0.0 .OR. pwy3 .NE. INT(pwy3))) & -& THEN - pwr3_d = 0.0_4 - ELSE - pwr3_d = pwy3*pwx3**(pwy3-1)*pwx3_d - END IF - pwr3 = pwx3**pwy3 - pr_imd_d = pwr3_d - ct*ht_d - ht*ct_d - pr_imd = pwr3 - ht*ct - ELSE - pr_imd_d = prr_d - pr_imd = prr - END IF - IF (1.e-6_sp .LT. ht + pr_imd/ct) THEN - ht_imd_d = ht_d + (pr_imd_d-pr_imd*ct_d/ct)/ct - ht_imd = ht + pr_imd/ct - ELSE - ht_imd = 1.e-6_sp - ht_imd_d = 0.0_4 - END IF - pwx1_d = ct*ht_imd_d + ht_imd*ct_d - pwx1 = ht_imd*ct - pwy1 = -nm1 - IF (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1))) & -& THEN - pwr1_d = 0.0_4 - ELSE - pwr1_d = pwy1*pwx1**(pwy1-1)*pwx1_d - END IF - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - IF (ct .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2))) THEN - pwr2_d = 0.0_4 - ELSE - pwr2_d = pwy2*ct**(pwy2-1)*ct_d - END IF - pwr2 = ct**pwy2 - pwx3_d = pwr1_d + pwr2_d - pwx3 = pwr1 + pwr2 - pwy3 = -d1pnm1 - IF (pwx3 .LE. 0.0 .AND. (pwy3 .EQ. 0.0 .OR. pwy3 .NE. INT(pwy3))) & -& THEN - pwr3_d = 0.0_4 - ELSE - pwr3_d = pwy3*pwx3**(pwy3-1)*pwx3_d - END IF - pwr3 = pwx3**pwy3 - ht_d = (pwr3_d-pwr3*ct_d/ct)/ct - ht = pwr3/ct - IF (0._sp .LT. prd + l) THEN - qd_d = prd_d + l_d - qd = prd + l - ELSE - qd = 0._sp - qd_d = 0.0_4 - END IF - qr_d = ct*(ht_imd_d-ht_d) + (ht_imd-ht)*ct_d - qr = (ht_imd-ht)*ct - q_d = qr_d + qd_d - q = qr + qd - END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG_D - -! Differentiation of gr_production_transfer_mlp_alg in reverse (adjoint) mode (with options fixinterface noISIZE context OpenMP) -!: -! gradient of useful results: q kexc hp ht en f_q cp pn ct -! with respect to varying inputs: kexc hp ht en f_q cp pn ct - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG_B(f_q, f_q_b, pn, pn_b, en, & -& en_b, cp, cp_b, beta, ct, ct_b, kexc, kexc_b, n, prcp, hp, hp_b, ht& -& , ht_b, q, q_b, pr, perc, l, prr, prd, qr, qd) - IMPLICIT NONE -! fixed NN output size - REAL(sp), DIMENSION(4), INTENT(IN) :: f_q - REAL(sp), DIMENSION(4) :: f_q_b - REAL(sp), INTENT(IN) :: pn, en, cp, beta, ct, kexc, n, prcp - REAL(sp) :: pn_b, en_b, cp_b, ct_b, kexc_b - REAL(sp), INTENT(INOUT) :: hp, ht, q - REAL(sp), INTENT(INOUT) :: hp_b, ht_b, q_b - REAL(sp) :: pr, perc, l, prr, prd, qr, qd - REAL(sp) :: inv_cp, ps, es, hp_imd, pr_imd, ht_imd, nm1, d1pnm1 - REAL(sp) :: inv_cp_b, ps_b, es_b, hp_imd_b, pr_imd_b, ht_imd_b - INTRINSIC TANH - INTRINSIC MAX - REAL(sp) :: pwx1 - REAL(sp) :: pwx1_b - REAL(sp) :: pwr1 - REAL(sp) :: pwr1_b - REAL(sp) :: pwy1 - REAL(sp) :: pwy2 - REAL(sp) :: pwr2 - REAL(sp) :: pwr2_b - REAL(sp) :: pwx3 - REAL(sp) :: pwx3_b - REAL(sp) :: pwy3 - REAL(sp) :: pwr3 - REAL(sp) :: pwr3_b - REAL(sp) :: temp - REAL(sp) :: temp0 - REAL(sp) :: temp_b - REAL(sp) :: temp1 - REAL(sp) :: temp2 - REAL(sp) :: temp3 - REAL(sp) :: temp_b0 - REAL(sp) :: temp_b1 - REAL(sp) :: temp4 - REAL(sp) :: temp_b2 - REAL(sp) :: temp_b3 - REAL(sp) :: temp_b4 - REAL(sp) :: temp_b5 - INTEGER :: branch - REAL(sp) :: perc_b - REAL(sp) :: pr_b - REAL(sp) :: qr_b - REAL(sp) :: prd_b - REAL(sp) :: l_b - REAL(sp) :: qd_b - REAL(sp) :: prr_b -!% Production - inv_cp = 1._sp/cp - pr = 0._sp - ps = cp*(1._sp-hp*hp)*TANH(pn*inv_cp)/(1._sp+hp*TANH(pn*inv_cp)) - CALL PUSHREAL4(ps) - ps = (1._sp+f_q(1))*ps - es = hp*cp*(2._sp-hp)*TANH(en*inv_cp)/(1._sp+(1._sp-hp)*TANH(en*& -& inv_cp)) - CALL PUSHREAL4(es) - es = (1._sp+f_q(2))*es - hp_imd = hp + (ps-es)*inv_cp - IF (pn .GT. 0) THEN - pr = pn - (hp_imd-hp)*cp - CALL PUSHCONTROL1B(0) - ELSE - CALL PUSHCONTROL1B(1) - END IF - pwx1 = 1._sp + (hp_imd/beta)**4 - pwr1 = pwx1**(-0.25_sp) - perc = hp_imd*cp*(1._sp-pwr1) -!% Transfer - nm1 = n - 1._sp - d1pnm1 = 1._sp/nm1 - l = (1._sp+f_q(4))*kexc*ht**3.5_sp - prr = 0.9_sp*(1._sp-f_q(3)**2)*(pr+perc) + l - prd = (0.1_sp+0.9_sp*f_q(3)**2)*(pr+perc) - IF (prcp .LT. 0._sp) THEN - pwx1 = ht*ct - pwy1 = -nm1 - CALL PUSHREAL4(pwr1) - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - pwr2 = ct**pwy2 - pwx3 = pwr1 - pwr2 - pwy3 = -d1pnm1 - pwr3 = pwx3**pwy3 - pr_imd = pwr3 - ht*ct - CALL PUSHCONTROL1B(1) - ELSE - pr_imd = prr - CALL PUSHCONTROL1B(0) - END IF - IF (1.e-6_sp .LT. ht + pr_imd/ct) THEN - ht_imd = ht + pr_imd/ct - CALL PUSHCONTROL1B(0) - ELSE - ht_imd = 1.e-6_sp - CALL PUSHCONTROL1B(1) - END IF - CALL PUSHREAL4(pwx1) - pwx1 = ht_imd*ct - CALL PUSHREAL4(pwy1) - pwy1 = -nm1 - CALL PUSHREAL4(pwr1) - pwr1 = pwx1**pwy1 - CALL PUSHREAL4(pwy2) - pwy2 = -nm1 - pwr2 = ct**pwy2 - CALL PUSHREAL4(pwx3) - pwx3 = pwr1 + pwr2 - CALL PUSHREAL4(pwy3) - pwy3 = -d1pnm1 - pwr3 = pwx3**pwy3 - CALL PUSHREAL4(ht) - ht = pwr3/ct - IF (0._sp .LT. prd + l) THEN - CALL PUSHCONTROL1B(0) - ELSE - CALL PUSHCONTROL1B(1) - END IF - pwx1 = ht_imd*ct - nm1 = n - 1._sp - pwy1 = -nm1 - pwy2 = -nm1 - d1pnm1 = 1._sp/nm1 - pwy3 = -d1pnm1 - inv_cp = 1._sp/cp - qr_b = q_b - qd_b = q_b - ht_imd_b = ct*qr_b -!$OMP ATOMIC update - ht_b = ht_b - ct*qr_b -!$OMP ATOMIC update - ct_b = ct_b + (ht_imd-ht)*qr_b - CALL POPCONTROL1B(branch) - IF (branch .EQ. 0) THEN - prd_b = qd_b - l_b = qd_b - ELSE - l_b = 0.0_4 - prd_b = 0.0_4 - END IF - CALL POPREAL4(ht) - pwr3_b = ht_b/ct - IF (pwx3 .LE. 0.0 .AND. (pwy3 .EQ. 0.0 .OR. pwy3 .NE. INT(pwy3))) & -& THEN - pwx3_b = 0.0_4 - ELSE - pwx3_b = pwy3*pwx3**(pwy3-1)*pwr3_b - END IF - CALL POPREAL4(pwy3) - CALL POPREAL4(pwx3) - pwr1_b = pwx3_b - pwr2_b = pwx3_b - IF (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1))) & -& THEN - pwx1_b = 0.0_4 - ELSE - pwx1_b = pwy1*pwx1**(pwy1-1)*pwr1_b - END IF - IF (ct .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2))) THEN -!$OMP ATOMIC update - ct_b = ct_b + ht_imd*pwx1_b - pwr3*ht_b/ct**2 - ELSE -!$OMP ATOMIC update - ct_b = ct_b + pwy2*ct**(pwy2-1)*pwr2_b - pwr3*ht_b/ct**2 + ht_imd*& -& pwx1_b - END IF - CALL POPREAL4(pwy2) - CALL POPREAL4(pwr1) - CALL POPREAL4(pwy1) - CALL POPREAL4(pwx1) - ht_imd_b = ht_imd_b + ct*pwx1_b - CALL POPCONTROL1B(branch) - IF (branch .EQ. 0) THEN - ht_b = ht_imd_b - pr_imd_b = ht_imd_b/ct -!$OMP ATOMIC update - ct_b = ct_b - pr_imd*ht_imd_b/ct**2 - ELSE - ht_b = 0.0_4 - pr_imd_b = 0.0_4 - END IF - CALL POPCONTROL1B(branch) - IF (branch .EQ. 0) THEN - prr_b = pr_imd_b - ELSE - pwr3_b = pr_imd_b - IF (pwx3 .LE. 0.0 .AND. (pwy3 .EQ. 0.0 .OR. pwy3 .NE. INT(pwy3))) & -& THEN - pwx3_b = 0.0_4 - ELSE - pwx3_b = pwy3*pwx3**(pwy3-1)*pwr3_b - END IF - pwr1_b = pwx3_b - pwr2_b = -pwx3_b - CALL POPREAL4(pwr1) - IF (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1))) & -& THEN - pwx1_b = 0.0_4 - ELSE - pwx1_b = pwy1*pwx1**(pwy1-1)*pwr1_b - END IF -!$OMP ATOMIC update - ht_b = ht_b + ct*pwx1_b - ct*pr_imd_b - IF (ct .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2))) & -& THEN -!$OMP ATOMIC update - ct_b = ct_b + ht*pwx1_b - ht*pr_imd_b - ELSE -!$OMP ATOMIC update - ct_b = ct_b + pwy2*ct**(pwy2-1)*pwr2_b - ht*pr_imd_b + ht*pwx1_b - END IF - pwx1 = 1._sp + (hp_imd/beta)**4 - prr_b = 0.0_4 - END IF -!$OMP ATOMIC update - f_q_b(3) = f_q_b(3) + 2*f_q(3)*0.9_sp*(pr+perc)*prd_b - 2*f_q(3)*(pr& -& +perc)*0.9_sp*prr_b - temp_b5 = (0.9_sp*f_q(3)**2+0.1_sp)*prd_b - pr_b = temp_b5 - perc_b = temp_b5 - temp_b5 = (1._sp-f_q(3)**2)*0.9_sp*prr_b - l_b = l_b + prr_b - pr_b = pr_b + temp_b5 - perc_b = perc_b + temp_b5 - inv_cp*hp_b - temp_b5 = ht**3.5_sp*l_b -!$OMP ATOMIC update - ht_b = ht_b + 3.5_sp*ht**2.5*(f_q(4)+1._sp)*kexc*l_b -!$OMP ATOMIC update - f_q_b(4) = f_q_b(4) + kexc*temp_b5 -!$OMP ATOMIC update - kexc_b = kexc_b + (f_q(4)+1._sp)*temp_b5 - inv_cp_b = -(perc*hp_b) -!$OMP ATOMIC update - cp_b = cp_b + hp_imd*(1._sp-pwr1)*perc_b - pwr1_b = -(hp_imd*cp*perc_b) - pwx1_b = -(0.25_sp*pwx1**(-1.25)*pwr1_b) - hp_imd_b = hp_b + cp*(1._sp-pwr1)*perc_b + 4*hp_imd**3*pwx1_b/beta**& -& 4 - CALL POPCONTROL1B(branch) - IF (branch .EQ. 0) THEN -!$OMP ATOMIC update - pn_b = pn_b + pr_b - hp_imd_b = hp_imd_b - cp*pr_b - hp_b = cp*pr_b -!$OMP ATOMIC update - cp_b = cp_b - (hp_imd-hp)*pr_b - ELSE - hp_b = 0.0_4 - END IF - es_b = -(inv_cp*hp_imd_b) - inv_cp_b = inv_cp_b + (ps-es)*hp_imd_b - CALL POPREAL4(es) -!$OMP ATOMIC update - f_q_b(2) = f_q_b(2) + es*es_b - es_b = (f_q(2)+1._sp)*es_b - temp4 = TANH(en*inv_cp) - temp3 = (-hp+1._sp)*temp4 + 1._sp - temp1 = TANH(en*inv_cp) - temp0 = hp*cp*(-hp+2._sp) - temp_b3 = es_b/temp3 - temp_b = (2._sp-hp)*temp1*temp_b3 - temp_b0 = -(temp0*temp1*temp_b3/temp3) -!$OMP ATOMIC update - hp_b = hp_b + hp_imd_b + cp*temp_b - hp*cp*temp1*temp_b3 - temp4*& -& temp_b0 - ps_b = inv_cp*hp_imd_b - temp_b4 = (1.0-TANH(en*inv_cp)**2)*temp0*temp_b3 - temp_b5 = (1.0-TANH(en*inv_cp)**2)*(1._sp-hp)*temp_b0 -!$OMP ATOMIC update - en_b = en_b + inv_cp*temp_b5 + inv_cp*temp_b4 -!$OMP ATOMIC update - cp_b = cp_b + hp*temp_b - CALL POPREAL4(ps) -!$OMP ATOMIC update - f_q_b(1) = f_q_b(1) + ps*ps_b - ps_b = (f_q(1)+1._sp)*ps_b - temp = TANH(pn*inv_cp) - temp0 = hp*temp + 1._sp - temp1 = TANH(pn*inv_cp) - temp2 = cp*(-(hp*hp)+1._sp) - temp_b = ps_b/temp0 - temp_b0 = (1.0-TANH(pn*inv_cp)**2)*temp2*temp_b - temp_b1 = -(temp2*temp1*temp_b/temp0) + be_b = be_b + EXP(ar)*qe_b + ar_b = EXP(ar)*be*qe_b + ELSE IF (branch .EQ. 1) THEN + ar = he_star/be + temp = EXP(ar) + he_star_b = he_star_b + qe_b !$OMP ATOMIC update - hp_b = hp_b + temp*temp_b1 - 2*hp*cp*temp1*temp_b - temp_b2 = (1.0-TANH(pn*inv_cp)**2)*hp*temp_b1 - inv_cp_b = inv_cp_b + en*temp_b5 + en*temp_b4 + pn*temp_b2 + pn*& -& temp_b0 + be_b = be_b + qe_b/temp + ar_b = -(EXP(ar)*be*qe_b/temp**2) + ELSE !$OMP ATOMIC update - cp_b = cp_b + (1._sp-hp**2)*temp1*temp_b - inv_cp_b/cp**2 + be_b = be_b + LOG(arg1)*qe_b + arg1_b = be*qe_b/arg1 + ar = he_star/be + ar_b = EXP(ar)*arg1_b + END IF + he_star_b = he_star_b + ar_b/be !$OMP ATOMIC update - pn_b = pn_b + inv_cp*temp_b2 + inv_cp*temp_b0 - END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG_B + be_b = be_b - he_star*ar_b/be**2 + he_b = he_star_b + pre_b = he_star_b + END SUBROUTINE GR_EXPONENTIAL_TRANSFER_B - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG(f_q, pn, en, cp, beta, ct, & -& kexc, n, prcp, hp, ht, q, pr, perc, l, prr, prd, qr, qd) + SUBROUTINE GR_EXPONENTIAL_TRANSFER(pre, be, he, qe) IMPLICIT NONE -! fixed NN output size - REAL(sp), DIMENSION(4), INTENT(IN) :: f_q - REAL(sp), INTENT(IN) :: pn, en, cp, beta, ct, kexc, n, prcp - REAL(sp), INTENT(INOUT) :: hp, ht, q - REAL(sp), INTENT(OUT) :: pr, perc, l, prr, prd, qr, qd - REAL(sp) :: inv_cp, ps, es, hp_imd, pr_imd, ht_imd, nm1, d1pnm1 - INTRINSIC TANH - INTRINSIC MAX - REAL(sp) :: pwx1 - REAL(sp) :: pwr1 - REAL(sp) :: pwy1 - REAL(sp) :: pwy2 - REAL(sp) :: pwr2 - REAL(sp) :: pwx3 - REAL(sp) :: pwy3 - REAL(sp) :: pwr3 -!% Production - inv_cp = 1._sp/cp - pr = 0._sp - ps = cp*(1._sp-hp*hp)*TANH(pn*inv_cp)/(1._sp+hp*TANH(pn*inv_cp)) - ps = (1._sp+f_q(1))*ps - es = hp*cp*(2._sp-hp)*TANH(en*inv_cp)/(1._sp+(1._sp-hp)*TANH(en*& -& inv_cp)) - es = (1._sp+f_q(2))*es - hp_imd = hp + (ps-es)*inv_cp - IF (pn .GT. 0) pr = pn - (hp_imd-hp)*cp - pwx1 = 1._sp + (hp_imd/beta)**4 - pwr1 = pwx1**(-0.25_sp) - perc = hp_imd*cp*(1._sp-pwr1) - hp = hp_imd - perc*inv_cp -!% Transfer - nm1 = n - 1._sp - d1pnm1 = 1._sp/nm1 - l = (1._sp+f_q(4))*kexc*ht**3.5_sp - prr = 0.9_sp*(1._sp-f_q(3)**2)*(pr+perc) + l - prd = (0.1_sp+0.9_sp*f_q(3)**2)*(pr+perc) - IF (prcp .LT. 0._sp) THEN - pwx1 = ht*ct - pwy1 = -nm1 - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - pwr2 = ct**pwy2 - pwx3 = pwr1 - pwr2 - pwy3 = -d1pnm1 - pwr3 = pwx3**pwy3 - pr_imd = pwr3 - ht*ct - ELSE - pr_imd = prr - END IF - IF (1.e-6_sp .LT. ht + pr_imd/ct) THEN - ht_imd = ht + pr_imd/ct - ELSE - ht_imd = 1.e-6_sp - END IF - pwx1 = ht_imd*ct - pwy1 = -nm1 - pwr1 = pwx1**pwy1 - pwy2 = -nm1 - pwr2 = ct**pwy2 - pwx3 = pwr1 + pwr2 - pwy3 = -d1pnm1 - pwr3 = pwx3**pwy3 - ht = pwr3/ct - IF (0._sp .LT. prd + l) THEN - qd = prd + l + REAL(sp), INTENT(IN) :: pre, be + REAL(sp), INTENT(INOUT) :: he + REAL(sp), INTENT(OUT) :: qe + REAL(sp) :: he_star, ar + INTRINSIC EXP + INTRINSIC LOG + REAL(sp) :: arg1 + he_star = he + pre + ar = he_star/be + IF (ar .LT. -7._sp) THEN + qe = be*EXP(ar) + ELSE IF (ar .GT. 7._sp) THEN + qe = he_star + be/EXP(ar) ELSE - qd = 0._sp + arg1 = EXP(ar) + 1._sp + qe = be*LOG(arg1) END IF - qr = (ht_imd-ht)*ct - q = qr + qd - END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ALG + he = he_star - qe + END SUBROUTINE GR_EXPONENTIAL_TRANSFER ! Differentiation of gr_production_transfer_ode in forward (tangent) mode (with options fixinterface noISIZE context OpenMP): ! variations of useful results: q hp ht @@ -14547,14 +14018,14 @@ END SUBROUTINE GR_PRODUCTION_TRANSFER_ODE ! Differentiation of gr_production_transfer_mlp_ode in forward (tangent) mode (with options fixinterface noISIZE context OpenMP) !: ! variations of useful results: q hp ht -! with respect to varying inputs: kexc hp ht en f_q cp pn ct - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(f_q, f_q_d, pn, pn_d, en, & +! with respect to varying inputs: kexc hp ht en fq cp pn ct + SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(fq, fq_d, pn, pn_d, en, & & en_d, cp, cp_d, ct, ct_d, kexc, kexc_d, hp, hp_d, ht, ht_d, q, q_d, & & l) IMPLICIT NONE ! fixed NN output size - REAL(sp), DIMENSION(5), INTENT(IN) :: f_q - REAL(sp), DIMENSION(5), INTENT(IN) :: f_q_d + REAL(sp), DIMENSION(5), INTENT(IN) :: fq + REAL(sp), DIMENSION(5), INTENT(IN) :: fq_d REAL(sp), INTENT(IN) :: pn, en, cp, ct, kexc REAL(sp), INTENT(IN) :: pn_d, en_d, cp_d, ct_d, kexc_d REAL(sp), INTENT(INOUT) :: hp, ht, q @@ -14573,11 +14044,11 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(f_q, f_q_d, pn, pn_d, en, & ! dt = 1._sp/real(n_subtimesteps, sp) dt = 1._sp !do i = 1, n_subtimesteps - temp = (f_q(2)+1._sp)*(-hp+2._sp) - fhp_d = (1._sp-hp**2)*(pn*f_q_d(1)+(f_q(1)+1._sp)*pn_d) - (f_q(1)+& -& 1._sp)*pn*2*hp*hp_d - en*hp*((2._sp-hp)*f_q_d(2)-(f_q(2)+1._sp)*& -& hp_d) - temp*(hp*en_d+en*hp_d) - fhp = (f_q(1)+1._sp)*pn*(1._sp-hp*hp) - temp*(en*hp) + temp = (fq(2)+1._sp)*(-hp+2._sp) + fhp_d = (1._sp-hp**2)*(pn*fq_d(1)+(fq(1)+1._sp)*pn_d) - (fq(1)+1._sp& +& )*pn*2*hp*hp_d - en*hp*((2._sp-hp)*fq_d(2)-(fq(2)+1._sp)*hp_d) - & +& temp*(hp*en_d+en*hp_d) + fhp = (fq(1)+1._sp)*pn*(1._sp-hp*hp) - temp*(en*hp) hp_d = hp_d + dt*(inv_cp*fhp_d+fhp*inv_cp_d) hp = hp + dt*fhp*inv_cp IF (hp .LE. 0._sp) THEN @@ -14588,16 +14059,16 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(f_q, f_q_d, pn, pn_d, en, & hp = 1._sp - 1.e-6_sp hp_d = 0.0_4 END IF - temp = (-(f_q(3)*f_q(3))+1._sp)*(hp*hp) + temp = (-(fq(3)*fq(3))+1._sp)*(hp*hp) temp0 = ht**3.5_sp temp1 = ht**5 - fht_d = 0.9_sp*((f_q(1)+1._sp)*pn*((1._sp-f_q(3)**2)*2*hp*hp_d-hp**2& -& *2*f_q(3)*f_q_d(3))+temp*(pn*f_q_d(1)+(f_q(1)+1._sp)*pn_d)) + & -& temp0*(kexc*f_q_d(4)+(f_q(4)+1._sp)*kexc_d) + ((f_q(4)+1._sp)*kexc& -& *3.5_sp*ht**2.5-(f_q(5)+1._sp)*ct*5*ht**4)*ht_d - temp1*(ct*f_q_d(& -& 5)+(f_q(5)+1._sp)*ct_d) - fht = 0.9_sp*(temp*((f_q(1)+1._sp)*pn)) + (f_q(4)+1._sp)*kexc*temp0 & -& - (f_q(5)+1._sp)*ct*temp1 + fht_d = 0.9_sp*((fq(1)+1._sp)*pn*((1._sp-fq(3)**2)*2*hp*hp_d-hp**2*2& +& *fq(3)*fq_d(3))+temp*(pn*fq_d(1)+(fq(1)+1._sp)*pn_d)) + temp0*(& +& kexc*fq_d(4)+(fq(4)+1._sp)*kexc_d) + ((fq(4)+1._sp)*kexc*3.5_sp*ht& +& **2.5-(fq(5)+1._sp)*ct*5*ht**4)*ht_d - temp1*(ct*fq_d(5)+(fq(5)+& +& 1._sp)*ct_d) + fht = 0.9_sp*(temp*((fq(1)+1._sp)*pn)) + (fq(4)+1._sp)*kexc*temp0 - & +& (fq(5)+1._sp)*ct*temp1 ht_d = ht_d + dt*(fht_d-fht*ct_d/ct)/ct ht = ht + dt*fht/ct IF (ht .LE. 0._sp) THEN @@ -14610,29 +14081,29 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D(f_q, f_q_d, pn, pn_d, en, & END IF !end do temp1 = ht**3.5_sp - l_d = temp1*(kexc*f_q_d(4)+(f_q(4)+1._sp)*kexc_d) + (f_q(4)+1._sp)*& -& kexc*3.5_sp*ht**2.5*ht_d - l = (f_q(4)+1._sp)*kexc*temp1 + l_d = temp1*(kexc*fq_d(4)+(fq(4)+1._sp)*kexc_d) + (fq(4)+1._sp)*kexc& +& *3.5_sp*ht**2.5*ht_d + l = (fq(4)+1._sp)*kexc*temp1 temp1 = ht**5 - temp0 = (f_q(1)+1._sp)*pn*(hp*hp) - temp = 0.9_sp*(f_q(3)*f_q(3)) + 0.1_sp - q_d = temp1*(ct*f_q_d(5)+(f_q(5)+1._sp)*ct_d) + (f_q(5)+1._sp)*ct*5*& -& ht**4*ht_d + temp0*0.9_sp*2*f_q(3)*f_q_d(3) + temp*(hp**2*(pn*& -& f_q_d(1)+(f_q(1)+1._sp)*pn_d)+(f_q(1)+1._sp)*pn*2*hp*hp_d) + l_d - q = (f_q(5)+1._sp)*ct*temp1 + temp*temp0 + l + temp0 = (fq(1)+1._sp)*pn*(hp*hp) + temp = 0.9_sp*(fq(3)*fq(3)) + 0.1_sp + q_d = temp1*(ct*fq_d(5)+(fq(5)+1._sp)*ct_d) + (fq(5)+1._sp)*ct*5*ht& +& **4*ht_d + temp0*0.9_sp*2*fq(3)*fq_d(3) + temp*(hp**2*(pn*fq_d(1)+& +& (fq(1)+1._sp)*pn_d)+(fq(1)+1._sp)*pn*2*hp*hp_d) + l_d + q = (fq(5)+1._sp)*ct*temp1 + temp*temp0 + l END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_D ! Differentiation of gr_production_transfer_mlp_ode in reverse (adjoint) mode (with options fixinterface noISIZE context OpenMP) !: -! gradient of useful results: q kexc hp ht en f_q cp pn ct -! with respect to varying inputs: kexc hp ht en f_q cp pn ct - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & +! gradient of useful results: q kexc hp ht en fq cp pn ct +! with respect to varying inputs: kexc hp ht en fq cp pn ct + SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(fq, fq_b, pn, pn_b, en, & & en_b, cp, cp_b, ct, ct_b, kexc, kexc_b, hp, hp_b, ht, ht_b, q, q_b, & & l) IMPLICIT NONE ! fixed NN output size - REAL(sp), DIMENSION(5), INTENT(IN) :: f_q - REAL(sp), DIMENSION(5) :: f_q_b + REAL(sp), DIMENSION(5), INTENT(IN) :: fq + REAL(sp), DIMENSION(5) :: fq_b REAL(sp), INTENT(IN) :: pn, en, cp, ct, kexc REAL(sp) :: pn_b, en_b, cp_b, ct_b, kexc_b REAL(sp), INTENT(INOUT) :: hp, ht, q @@ -14652,8 +14123,8 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & ! dt = 1._sp/real(n_subtimesteps, sp) dt = 1._sp !do i = 1, n_subtimesteps - fhp = (1._sp+f_q(1))*pn*(1._sp-hp**2) - (1._sp+f_q(2))*en*hp*(2._sp-& -& hp) + fhp = (1._sp+fq(1))*pn*(1._sp-hp**2) - (1._sp+fq(2))*en*hp*(2._sp-hp& +& ) CALL PUSHREAL4(hp) hp = hp + dt*fhp*inv_cp IF (hp .LE. 0._sp) THEN @@ -14668,8 +14139,8 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & ELSE CALL PUSHCONTROL1B(1) END IF - fht = 0.9_sp*(1._sp-f_q(3)**2)*(1._sp+f_q(1))*pn*hp**2 + (1._sp+f_q(& -& 4))*kexc*ht**3.5_sp - (1._sp+f_q(5))*ct*ht**5 + fht = 0.9_sp*(1._sp-fq(3)**2)*(1._sp+fq(1))*pn*hp**2 + (1._sp+fq(4))& +& *kexc*ht**3.5_sp - (1._sp+fq(5))*ct*ht**5 CALL PUSHREAL4(ht) ht = ht + dt*fht/ct IF (ht .LE. 0._sp) THEN @@ -14687,27 +14158,27 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & l_b = q_b temp_b2 = ht**5*q_b !$OMP ATOMIC update - ht_b = ht_b + 5*ht**4*(f_q(5)+1._sp)*ct*q_b + 3.5_sp*ht**2.5*(f_q(4)& -& +1._sp)*kexc*l_b + ht_b = ht_b + 5*ht**4*(fq(5)+1._sp)*ct*q_b + 3.5_sp*ht**2.5*(fq(4)+& +& 1._sp)*kexc*l_b !$OMP ATOMIC update - f_q_b(3) = f_q_b(3) + 2*f_q(3)*0.9_sp*(f_q(1)+1._sp)*pn*hp**2*q_b - temp_b1 = (0.9_sp*f_q(3)**2+0.1_sp)*q_b + fq_b(3) = fq_b(3) + 2*fq(3)*0.9_sp*(fq(1)+1._sp)*pn*hp**2*q_b + temp_b1 = (0.9_sp*fq(3)**2+0.1_sp)*q_b temp_b0 = hp**2*temp_b1 !$OMP ATOMIC update - hp_b = hp_b + 2*hp*(f_q(1)+1._sp)*pn*temp_b1 + hp_b = hp_b + 2*hp*(fq(1)+1._sp)*pn*temp_b1 !$OMP ATOMIC update - f_q_b(1) = f_q_b(1) + pn*temp_b0 + fq_b(1) = fq_b(1) + pn*temp_b0 !$OMP ATOMIC update - pn_b = pn_b + (f_q(1)+1._sp)*temp_b0 + pn_b = pn_b + (fq(1)+1._sp)*temp_b0 !$OMP ATOMIC update - f_q_b(5) = f_q_b(5) + ct*temp_b2 + fq_b(5) = fq_b(5) + ct*temp_b2 !$OMP ATOMIC update - ct_b = ct_b + (f_q(5)+1._sp)*temp_b2 + ct_b = ct_b + (fq(5)+1._sp)*temp_b2 temp_b2 = ht**3.5_sp*l_b !$OMP ATOMIC update - f_q_b(4) = f_q_b(4) + kexc*temp_b2 + fq_b(4) = fq_b(4) + kexc*temp_b2 !$OMP ATOMIC update - kexc_b = kexc_b + (f_q(4)+1._sp)*temp_b2 + kexc_b = kexc_b + (fq(4)+1._sp)*temp_b2 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) ht_b = 0.0_4 CALL POPCONTROL1B(branch) @@ -14718,29 +14189,29 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & fht_b = temp_b2 !$OMP ATOMIC update ct_b = ct_b - fht*temp_b2/ct - temp_b1 = (f_q(1)+1._sp)*pn*0.9_sp*fht_b - temp_b0 = (1._sp-f_q(3)**2)*hp**2*0.9_sp*fht_b + temp_b1 = (fq(1)+1._sp)*pn*0.9_sp*fht_b + temp_b0 = (1._sp-fq(3)**2)*hp**2*0.9_sp*fht_b temp_b = ht**3.5_sp*fht_b !$OMP ATOMIC update - ht_b = ht_b + (3.5_sp*ht**2.5*(f_q(4)+1._sp)*kexc-5*ht**4*(f_q(5)+& + ht_b = ht_b + (3.5_sp*ht**2.5*(fq(4)+1._sp)*kexc-5*ht**4*(fq(5)+& & 1._sp)*ct)*fht_b temp_b2 = -(ht**5*fht_b) !$OMP ATOMIC update - f_q_b(5) = f_q_b(5) + ct*temp_b2 + fq_b(5) = fq_b(5) + ct*temp_b2 !$OMP ATOMIC update - ct_b = ct_b + (f_q(5)+1._sp)*temp_b2 + ct_b = ct_b + (fq(5)+1._sp)*temp_b2 !$OMP ATOMIC update - f_q_b(4) = f_q_b(4) + kexc*temp_b + fq_b(4) = fq_b(4) + kexc*temp_b !$OMP ATOMIC update - kexc_b = kexc_b + (f_q(4)+1._sp)*temp_b + kexc_b = kexc_b + (fq(4)+1._sp)*temp_b !$OMP ATOMIC update - f_q_b(1) = f_q_b(1) + pn*temp_b0 + fq_b(1) = fq_b(1) + pn*temp_b0 !$OMP ATOMIC update - pn_b = pn_b + (f_q(1)+1._sp)*temp_b0 + pn_b = pn_b + (fq(1)+1._sp)*temp_b0 !$OMP ATOMIC update - f_q_b(3) = f_q_b(3) - 2*f_q(3)*hp**2*temp_b1 + fq_b(3) = fq_b(3) - 2*fq(3)*hp**2*temp_b1 !$OMP ATOMIC update - hp_b = hp_b + 2*hp*(1._sp-f_q(3)**2)*temp_b1 + hp_b = hp_b + 2*hp*(1._sp-fq(3)**2)*temp_b1 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) hp_b = 0.0_4 CALL POPCONTROL1B(branch) @@ -14751,27 +14222,27 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B(f_q, f_q_b, pn, pn_b, en, & inv_cp_b = fhp*dt*hp_b temp_b = (1._sp-hp**2)*fhp_b temp_b0 = -(en*hp*fhp_b) - temp_b1 = -((f_q(2)+1._sp)*(2._sp-hp)*fhp_b) + temp_b1 = -((fq(2)+1._sp)*(2._sp-hp)*fhp_b) !$OMP ATOMIC update - hp_b = hp_b + en*temp_b1 - 2*hp*(f_q(1)+1._sp)*pn*fhp_b - (f_q(2)+& + hp_b = hp_b + en*temp_b1 - 2*hp*(fq(1)+1._sp)*pn*fhp_b - (fq(2)+& & 1._sp)*temp_b0 !$OMP ATOMIC update en_b = en_b + hp*temp_b1 !$OMP ATOMIC update - f_q_b(2) = f_q_b(2) + (2._sp-hp)*temp_b0 + fq_b(2) = fq_b(2) + (2._sp-hp)*temp_b0 !$OMP ATOMIC update - f_q_b(1) = f_q_b(1) + pn*temp_b + fq_b(1) = fq_b(1) + pn*temp_b !$OMP ATOMIC update - pn_b = pn_b + (f_q(1)+1._sp)*temp_b + pn_b = pn_b + (fq(1)+1._sp)*temp_b !$OMP ATOMIC update cp_b = cp_b - inv_cp_b/cp**2 END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE_B - SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE(f_q, pn, en, cp, ct, kexc, & -& hp, ht, q, l) + SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE(fq, pn, en, cp, ct, kexc, hp& +& , ht, q, l) IMPLICIT NONE ! fixed NN output size - REAL(sp), DIMENSION(5), INTENT(IN) :: f_q + REAL(sp), DIMENSION(5), INTENT(IN) :: fq REAL(sp), INTENT(IN) :: pn, en, cp, ct, kexc REAL(sp), INTENT(INOUT) :: hp, ht, q REAL(sp), INTENT(OUT) :: l @@ -14782,20 +14253,20 @@ SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE(f_q, pn, en, cp, ct, kexc, & ! dt = 1._sp/real(n_subtimesteps, sp) dt = 1._sp !do i = 1, n_subtimesteps - fhp = (1._sp+f_q(1))*pn*(1._sp-hp**2) - (1._sp+f_q(2))*en*hp*(2._sp-& -& hp) + fhp = (1._sp+fq(1))*pn*(1._sp-hp**2) - (1._sp+fq(2))*en*hp*(2._sp-hp& +& ) hp = hp + dt*fhp*inv_cp IF (hp .LE. 0._sp) hp = 1.e-6_sp IF (hp .GE. 1._sp) hp = 1._sp - 1.e-6_sp - fht = 0.9_sp*(1._sp-f_q(3)**2)*(1._sp+f_q(1))*pn*hp**2 + (1._sp+f_q(& -& 4))*kexc*ht**3.5_sp - (1._sp+f_q(5))*ct*ht**5 + fht = 0.9_sp*(1._sp-fq(3)**2)*(1._sp+fq(1))*pn*hp**2 + (1._sp+fq(4))& +& *kexc*ht**3.5_sp - (1._sp+fq(5))*ct*ht**5 ht = ht + dt*fht/ct IF (ht .LE. 0._sp) ht = 1.e-6_sp IF (ht .GE. 1._sp) ht = 1._sp - 1.e-6_sp !end do - l = (1._sp+f_q(4))*kexc*ht**3.5_sp - q = (1._sp+f_q(5))*ct*ht**5 + (0.1_sp+0.9_sp*f_q(3)**2)*(1._sp+f_q(1& -& ))*pn*hp**2 + l + l = (1._sp+fq(4))*kexc*ht**3.5_sp + q = (1._sp+fq(5))*ct*ht**5 + (0.1_sp+0.9_sp*fq(3)**2)*(1._sp+fq(1))*& +& pn*hp**2 + l END SUBROUTINE GR_PRODUCTION_TRANSFER_MLP_ODE ! Differentiation of gr4_time_step in forward (tangent) mode (with options fixinterface noISIZE context OpenMP): @@ -14855,11 +14326,11 @@ SUBROUTINE GR4_TIME_STEP_D(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION_D(ac_prcp(k), ac_prcp_d(k), ac_pet(k), & & ac_ci(k), ac_ci_d(k), ac_hi(k), ac_hi_d(k)& & , pn, pn_d, en, en_d) - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k& -& ), beta, ac_hp(k), ac_hp_d(k), pr, pr_d, perc& -& , perc_d) - CALL GR_EXCHANGE_D(ac_kexc(k), ac_kexc_d(k), ac_ht(k), & -& ac_ht_d(k), l, l_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_cp(k), ac_cp_d(k), beta, ac_hp(k& +& ), ac_hp_d(k), pr, pr_d, perc, perc_d) + CALL GR_EXCHANGE_D(0._sp, 0.0_4, ac_kexc(k), ac_kexc_d(k), & +& ac_ht(k), ac_ht_d(k), l, l_d) ELSE pr = 0._sp perc = 0._sp @@ -14925,6 +14396,9 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, returns, & REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd REAL(sp) :: pn_b, en_b, pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 + REAL(sp) :: dummydiff_b1 INTEGER :: branch INTEGER :: chunk_start INTEGER :: chunk_end @@ -14955,9 +14429,9 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) - CALL GR_EXCHANGE(ac_kexc(k), ac_ht(k), l) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) + CALL GR_EXCHANGE(0._sp, ac_kexc(k), ac_ht(k), l) CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) @@ -15032,12 +14506,15 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, returns, & l_b = l_b + prr_b CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN - CALL GR_EXCHANGE_B(ac_kexc(k), ac_kexc_b(k), ac_ht(k), & -& ac_ht_b(k), l, l_b) + CALL GR_EXCHANGE_B(0._sp, dummydiff_b1, ac_kexc(k), & +& ac_kexc_b(k), ac_ht(k), ac_ht_b(k), l, l_b) CALL POPREAL4(ac_hp(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k& -& ), beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc& -& , perc_b) + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k), & +& beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc, & +& perc_b) CALL POPREAL4(ac_hi(k)) CALL POPREAL4(pn) CALL POPREAL4(en) @@ -15091,9 +14568,9 @@ SUBROUTINE GR4_TIME_STEP(setup, mesh, input_data, options, returns, & IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) - CALL GR_EXCHANGE(ac_kexc(k), ac_ht(k), l) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) + CALL GR_EXCHANGE(0._sp, ac_kexc(k), ac_ht(k), l) ELSE pr = 0._sp perc = 0._sp @@ -15167,6 +14644,8 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_D(setup, mesh, input_data, options, & REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_d, pn_d, en_d INTEGER :: row, col, k, time_step_returns REAL(sp) :: beta, pr, perc, l, prr, prd, qr, qd + REAL(sp) :: pr_d, perc_d, l_d, prr_d, prd_d, qr_d, qd_d + INTRINSIC MAX REAL(sp) :: temp CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) @@ -15201,50 +14680,76 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_D(setup, mesh, input_data, options, & END DO END DO output_layer_d = 0.0_4 - input_layer_d = 0.0_4 ! Forward MLP without OPENMP DO col=1,mesh%ncol DO row=1,mesh%nrow IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - input_layer_d(1) = ac_hp_d(k) - input_layer(1) = ac_hp(k) - input_layer_d(2) = ac_ht_d(k) - input_layer(2) = ac_ht(k) - input_layer_d(3) = pn_d(k) - input_layer(3) = pn(k) - input_layer_d(4) = en_d(k) - input_layer(4) = en(k) - CALL FORWARD_MLP_D(weight_1, weight_1_d, bias_1, bias_1_d, & -& weight_2, weight_2_d, bias_2, bias_2_d, & -& input_layer, input_layer_d, output_layer(:, k), & -& output_layer_d(:, k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + input_layer_d(:) = (/ac_hp_d(k), ac_ht_d(k), pn_d(k), en_d(k& +& )/) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP_D(weight_1, weight_1_d, bias_1, bias_1_d, & +& weight_2, weight_2_d, bias_2, bias_2_d, & +& input_layer, input_layer_d, output_layer(:, k)& +& , output_layer_d(:, k)) + ELSE + output_layer_d(:, k) = 0.0_4 + output_layer(:, k) = 0._sp + END IF END IF END DO END DO ! Production and transfer with OPENMP !$OMP PARALLEL DO NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, & -!$OMP&returns, output_layer, ac_prcp, ac_cp, beta, ac_ct, ac_kexc, ac_hp& -!$OMP&, ac_ht, ac_qt, pn, en), SHARED(output_layer_d, ac_prcp_d, ac_cp_d& -!$OMP&, ac_ct_d, ac_kexc_d, ac_hp_d, ac_ht_d, ac_qt_d, pn_d, en_d), & -!$OMP&PRIVATE(row, col, k, time_step_returns, pr, perc, l, prr, prd, qr& -!$OMP&, qd), PRIVATE(temp), SCHEDULE(static) +!$OMP&returns, output_layer, ac_prcp, ac_pet, ac_cp, beta, ac_ct, & +!$OMP&ac_kexc, ac_hp, ac_ht, ac_qt, pn, en), SHARED(output_layer_d, & +!$OMP&ac_prcp_d, ac_cp_d, ac_ct_d, ac_kexc_d, ac_hp_d, ac_ht_d, ac_qt_d& +!$OMP&, pn_d, en_d), PRIVATE(row, col, k, time_step_returns, pr, perc, l& +!$OMP&, prr, prd, qr, qd), PRIVATE(pr_d, perc_d, l_d, prr_d, prd_d, qr_d& +!$OMP&, qd_d), PRIVATE(temp), SCHEDULE(static) DO col=1,mesh%ncol DO row=1,mesh%nrow IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - CALL GR_PRODUCTION_TRANSFER_MLP_ALG_D(output_layer(:, k), & -& output_layer_d(:, k), pn(k), & -& pn_d(k), en(k), en_d(k), ac_cp& -& (k), ac_cp_d(k), beta, ac_ct(k& -& ), ac_ct_d(k), ac_kexc(k), & -& ac_kexc_d(k), 5._sp, ac_prcp(k& -& ), ac_hp(k), ac_hp_d(k), ac_ht& -& (k), ac_ht_d(k), ac_qt(k), & -& ac_qt_d(k), pr, perc, l, prr, & -& prd, qr, qd) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL GR_PRODUCTION_D(output_layer(1, k), output_layer_d(1, k& +& ), output_layer(2, k), output_layer_d(2, k), & +& pn(k), pn_d(k), en(k), en_d(k), ac_cp(k), & +& ac_cp_d(k), beta, ac_hp(k), ac_hp_d(k), pr, & +& pr_d, perc, perc_d) + CALL GR_EXCHANGE_D(output_layer(4, k), output_layer_d(4, k)& +& , ac_kexc(k), ac_kexc_d(k), ac_ht(k), ac_ht_d(k& +& ), l, l_d) + ELSE + pr = 0._sp + perc = 0._sp + l = 0._sp + l_d = 0.0_4 + perc_d = 0.0_4 + pr_d = 0.0_4 + END IF + temp = -(output_layer(3, k)*output_layer(3, k)) + 1._sp + prr_d = 0.9_sp*(temp*(pr_d+perc_d)-(pr+perc)*2*output_layer(3& +& , k)*output_layer_d(3, k)) + l_d + prr = 0.9_sp*(temp*(pr+perc)) + l + temp = 0.9_sp*(output_layer(3, k)*output_layer(3, k)) + 0.1_sp + prd_d = (pr+perc)*0.9_sp*2*output_layer(3, k)*output_layer_d(3& +& , k) + temp*(pr_d+perc_d) + prd = temp*(pr+perc) + CALL GR_TRANSFER_D(5._sp, ac_prcp(k), prr, prr_d, ac_ct(k), & +& ac_ct_d(k), ac_ht(k), ac_ht_d(k), qr, qr_d) + IF (0._sp .LT. prd + l) THEN + qd_d = prd_d + l_d + qd = prd + l + ELSE + qd = 0._sp + qd_d = 0.0_4 + END IF + ac_qt_d(k) = qr_d + qd_d + ac_qt(k) = qr + qd ! Transform from mm/dt to m3/s temp = 1e-3_sp*mesh%dx(row, col)*mesh%dy(row, col) ac_qt_d(k) = temp*ac_qt_d(k)/setup%dt @@ -15306,9 +14811,12 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_B(setup, mesh, input_data, options, & REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_b, pn_b, en_b INTEGER :: row, col, k, time_step_returns REAL(sp) :: beta, pr, perc, l, prr, prd, qr, qd + REAL(sp) :: pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b + INTRINSIC MAX INTEGER :: branch INTEGER :: chunk_start INTEGER :: chunk_end + REAL(sp) :: temp_b CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& @@ -15347,28 +14855,28 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_B(setup, mesh, input_data, options, & DO row=1,mesh%nrow IF (mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0) THEN - CALL PUSHCONTROL1B(0) + CALL PUSHCONTROL2B(0) ELSE k = mesh%rowcol_to_ind_ac(row, col) - CALL PUSHREAL4(input_layer(1)) - input_layer(1) = ac_hp(k) - CALL PUSHREAL4(input_layer(2)) - input_layer(2) = ac_ht(k) - CALL PUSHREAL4(input_layer(3)) - input_layer(3) = pn(k) - CALL PUSHREAL4(input_layer(4)) - input_layer(4) = en(k) - CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & -& input_layer, output_layer(:, k)) - CALL PUSHCONTROL1B(1) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL PUSHREAL4ARRAY(input_layer, 4) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & +& input_layer, output_layer(:, k)) + CALL PUSHCONTROL2B(2) + ELSE + output_layer(:, k) = 0._sp + CALL PUSHCONTROL2B(1) + END IF END IF END DO END DO ! Production and transfer with OPENMP !$OMP PARALLEL NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, & -!$OMP&returns, output_layer, ac_prcp, ac_cp, beta, ac_ct, ac_kexc, ac_hp& -!$OMP&, ac_ht, ac_qt, pn, en), PRIVATE(row, col, k, time_step_returns, & -!$OMP&pr, perc, l, prr, prd, qr, qd), PRIVATE(chunk_start, chunk_end) +!$OMP&returns, output_layer, ac_prcp, ac_pet, ac_cp, beta, ac_ct, & +!$OMP&ac_kexc, ac_hp, ac_ht, ac_qt, pn, en), PRIVATE(row, col, k, & +!$OMP&time_step_returns, pr, perc, l, prr, prd, qr, qd), PRIVATE(& +!$OMP&chunk_start, chunk_end) CALL GETSTATICSCHEDULE(1, mesh%ncol, 1, chunk_start, chunk_end) DO col=chunk_start,chunk_end DO row=1,mesh%nrow @@ -15376,79 +14884,142 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP_B(setup, mesh, input_data, options, & & local_active_cell(row, col) .EQ. 0) THEN CALL PUSHCONTROL1B(0) ELSE + CALL PUSHINTEGER4(k) k = mesh%rowcol_to_ind_ac(row, col) - CALL PUSHREAL4(ac_qt(k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL PUSHREAL4(perc) + CALL PUSHREAL4(pr) + CALL PUSHREAL4(ac_hp(k)) + CALL GR_PRODUCTION(output_layer(1, k), output_layer(2, k), & +& pn(k), en(k), ac_cp(k), beta, ac_hp(k), pr, & +& perc) + CALL GR_EXCHANGE(output_layer(4, k), ac_kexc(k), ac_ht(k), l& +& ) + CALL PUSHCONTROL1B(1) + ELSE + CALL PUSHREAL4(pr) + pr = 0._sp + CALL PUSHREAL4(perc) + perc = 0._sp + l = 0._sp + CALL PUSHCONTROL1B(0) + END IF + CALL PUSHREAL4(prr) + prr = 0.9_sp*(1._sp-output_layer(3, k)**2)*(pr+perc) + l + prd = (0.1_sp+0.9_sp*output_layer(3, k)**2)*(pr+perc) CALL PUSHREAL4(ac_ht(k)) - CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION_TRANSFER_MLP_ALG(output_layer(:, k), pn(k)& -& , en(k), ac_cp(k), beta, ac_ct(k& -& ), ac_kexc(k), 5._sp, ac_prcp(k)& -& , ac_hp(k), ac_ht(k), ac_qt(k), & -& pr, perc, l, prr, prd, qr, qd) -! Transform from mm/dt to m3/s + CALL GR_TRANSFER(5._sp, ac_prcp(k), prr, ac_ct(k), ac_ht(k), & +& qr) + IF (0._sp .LT. prd + l) THEN + CALL PUSHCONTROL1B(0) + ELSE + CALL PUSHCONTROL1B(1) + END IF CALL PUSHCONTROL1B(1) END IF END DO END DO + CALL PUSHREAL4(pr) + CALL PUSHREAL4(perc) + CALL PUSHREAL4(prr) + CALL PUSHINTEGER4(k) !$OMP END PARALLEL output_layer_b = 0.0_4 en_b = 0.0_4 pn_b = 0.0_4 !$OMP PARALLEL NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, & -!$OMP&returns, output_layer, ac_prcp, ac_cp, beta, ac_ct, ac_kexc, ac_hp& -!$OMP&, ac_ht, ac_qt, pn, en), SHARED(output_layer_b, ac_prcp_b, ac_cp_b& -!$OMP&, ac_ct_b, ac_kexc_b, ac_hp_b, ac_ht_b, ac_qt_b, pn_b, en_b), & -!$OMP&PRIVATE(row, col, k, time_step_returns, pr, perc, l, prr, prd, qr& -!$OMP&, qd), PRIVATE(branch, chunk_end, chunk_start) +!$OMP&returns, output_layer, ac_prcp, ac_pet, ac_cp, beta, ac_ct, & +!$OMP&ac_kexc, ac_hp, ac_ht, ac_qt, pn, en), SHARED(output_layer_b, & +!$OMP&ac_prcp_b, ac_cp_b, ac_ct_b, ac_kexc_b, ac_hp_b, ac_ht_b, ac_qt_b& +!$OMP&, pn_b, en_b), PRIVATE(row, col, k, time_step_returns, pr, perc, l& +!$OMP&, prr, prd, qr, qd), PRIVATE(pr_b, perc_b, l_b, prr_b, prd_b, qr_b& +!$OMP&, qd_b), PRIVATE(branch, chunk_end, chunk_start), PRIVATE(temp_b) + CALL POPINTEGER4(k) + CALL POPREAL4(prr) + CALL POPREAL4(perc) + CALL POPREAL4(pr) + pr_b = 0.0_4 + perc_b = 0.0_4 + l_b = 0.0_4 + prr_b = 0.0_4 + prd_b = 0.0_4 + qr_b = 0.0_4 + qd_b = 0.0_4 CALL GETSTATICSCHEDULE(1, mesh%ncol, 1, chunk_start, chunk_end) DO col=chunk_end,chunk_start,-1 DO row=mesh%nrow,1,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN - k = mesh%rowcol_to_ind_ac(row, col) ac_qt_b(k) = mesh%dx(row, col)*1e-3_sp*mesh%dy(row, col)*& & ac_qt_b(k)/setup%dt - CALL POPREAL4(ac_hp(k)) - CALL POPREAL4(ac_ht(k)) - CALL POPREAL4(ac_qt(k)) - CALL GR_PRODUCTION_TRANSFER_MLP_ALG_B(output_layer(:, k), & -& output_layer_b(:, k), pn(k), & -& pn_b(k), en(k), en_b(k), ac_cp& -& (k), ac_cp_b(k), beta, ac_ct(k& -& ), ac_ct_b(k), ac_kexc(k), & -& ac_kexc_b(k), 5._sp, ac_prcp(k& -& ), ac_hp(k), ac_hp_b(k), ac_ht& -& (k), ac_ht_b(k), ac_qt(k), & -& ac_qt_b(k), pr, perc, l, prr, & -& prd, qr, qd) + qr_b = ac_qt_b(k) + qd_b = ac_qt_b(k) ac_qt_b(k) = 0.0_4 + CALL POPCONTROL1B(branch) + IF (branch .EQ. 0) THEN + prd_b = qd_b + l_b = qd_b + ELSE + l_b = 0.0_4 + prd_b = 0.0_4 + END IF + CALL POPREAL4(ac_ht(k)) + CALL GR_TRANSFER_B(5._sp, ac_prcp(k), prr, prr_b, ac_ct(k), & +& ac_ct_b(k), ac_ht(k), ac_ht_b(k), qr, qr_b) +!$OMP ATOMIC update + output_layer_b(3, k) = output_layer_b(3, k) + 2*output_layer(3& +& , k)*0.9_sp*(pr+perc)*prd_b - 2*output_layer(3, k)*(pr+perc)& +& *0.9_sp*prr_b + temp_b = (0.9_sp*output_layer(3, k)**2+0.1_sp)*prd_b + pr_b = temp_b + perc_b = temp_b + CALL POPREAL4(prr) + temp_b = (1._sp-output_layer(3, k)**2)*0.9_sp*prr_b + l_b = l_b + prr_b + pr_b = pr_b + temp_b + perc_b = perc_b + temp_b + CALL POPCONTROL1B(branch) + IF (branch .EQ. 0) THEN + CALL POPREAL4(perc) + CALL POPREAL4(pr) + ELSE + CALL GR_EXCHANGE_B(output_layer(4, k), output_layer_b(4, k)& +& , ac_kexc(k), ac_kexc_b(k), ac_ht(k), ac_ht_b(k& +& ), l, l_b) + CALL POPREAL4(ac_hp(k)) + CALL POPREAL4(pr) + CALL POPREAL4(perc) + CALL GR_PRODUCTION_B(output_layer(1, k), output_layer_b(1, k& +& ), output_layer(2, k), output_layer_b(2, k), & +& pn(k), pn_b(k), en(k), en_b(k), ac_cp(k), & +& ac_cp_b(k), beta, ac_hp(k), ac_hp_b(k), pr, & +& pr_b, perc, perc_b) + END IF + CALL POPINTEGER4(k) END IF END DO END DO !$OMP END PARALLEL - input_layer_b = 0.0_4 DO col=mesh%ncol,1,-1 DO row=mesh%nrow,1,-1 - CALL POPCONTROL1B(branch) + CALL POPCONTROL2B(branch) IF (branch .NE. 0) THEN - k = mesh%rowcol_to_ind_ac(row, col) - CALL FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & -& weight_2, weight_2_b, bias_2, bias_2_b, & -& input_layer, input_layer_b, output_layer(:, k), & -& output_layer_b(:, k)) - output_layer_b(:, k) = 0.0_4 - CALL POPREAL4(input_layer(4)) - en_b(k) = en_b(k) + input_layer_b(4) - input_layer_b(4) = 0.0_4 - CALL POPREAL4(input_layer(3)) - pn_b(k) = pn_b(k) + input_layer_b(3) - input_layer_b(3) = 0.0_4 - CALL POPREAL4(input_layer(2)) - ac_ht_b(k) = ac_ht_b(k) + input_layer_b(2) - input_layer_b(2) = 0.0_4 - CALL POPREAL4(input_layer(1)) - ac_hp_b(k) = ac_hp_b(k) + input_layer_b(1) - input_layer_b(1) = 0.0_4 + IF (branch .EQ. 1) THEN + k = mesh%rowcol_to_ind_ac(row, col) + output_layer_b(:, k) = 0.0_4 + ELSE + k = mesh%rowcol_to_ind_ac(row, col) + CALL FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & +& weight_2, weight_2_b, bias_2, bias_2_b, & +& input_layer, input_layer_b, output_layer(:, k)& +& , output_layer_b(:, k)) + output_layer_b(:, k) = 0.0_4 + CALL POPREAL4ARRAY(input_layer, 4) + ac_hp_b(k) = ac_hp_b(k) + input_layer_b(1) + ac_ht_b(k) = ac_ht_b(k) + input_layer_b(2) + pn_b(k) = pn_b(k) + input_layer_b(3) + en_b(k) = en_b(k) + input_layer_b(4) + END IF END IF END DO END DO @@ -15512,6 +15083,7 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP(setup, mesh, input_data, options, & REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet, pn, en INTEGER :: row, col, k, time_step_returns REAL(sp) :: beta, pr, perc, l, prr, prd, qr, qd + INTRINSIC MAX CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& & , 'prcp', ac_prcp) CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step& @@ -15544,30 +15116,47 @@ SUBROUTINE GR4_MLP_ALG_TIME_STEP(setup, mesh, input_data, options, & IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - input_layer(1) = ac_hp(k) - input_layer(2) = ac_ht(k) - input_layer(3) = pn(k) - input_layer(4) = en(k) - CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & -& input_layer, output_layer(:, k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & +& input_layer, output_layer(:, k)) + ELSE + output_layer(:, k) = 0._sp + END IF END IF END DO END DO ! Production and transfer with OPENMP !$OMP PARALLEL DO NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, & -!$OMP&returns, output_layer, ac_prcp, ac_cp, beta, ac_ct, ac_kexc, ac_hp& -!$OMP&, ac_ht, ac_qt, pn, en), PRIVATE(row, col, k, time_step_returns, & -!$OMP&pr, perc, l, prr, prd, qr, qd), SCHEDULE(static) +!$OMP&returns, output_layer, ac_prcp, ac_pet, ac_cp, beta, ac_ct, & +!$OMP&ac_kexc, ac_hp, ac_ht, ac_qt, pn, en), PRIVATE(row, col, k, & +!$OMP&time_step_returns, pr, perc, l, prr, prd, qr, qd), SCHEDULE(static) DO col=1,mesh%ncol DO row=1,mesh%nrow IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - CALL GR_PRODUCTION_TRANSFER_MLP_ALG(output_layer(:, k), pn(k)& -& , en(k), ac_cp(k), beta, ac_ct(k& -& ), ac_kexc(k), 5._sp, ac_prcp(k)& -& , ac_hp(k), ac_ht(k), ac_qt(k), & -& pr, perc, l, prr, prd, qr, qd) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL GR_PRODUCTION(output_layer(1, k), output_layer(2, k), & +& pn(k), en(k), ac_cp(k), beta, ac_hp(k), pr, & +& perc) + CALL GR_EXCHANGE(output_layer(4, k), ac_kexc(k), ac_ht(k), l& +& ) + ELSE + pr = 0._sp + perc = 0._sp + l = 0._sp + END IF + prr = 0.9_sp*(1._sp-output_layer(3, k)**2)*(pr+perc) + l + prd = (0.1_sp+0.9_sp*output_layer(3, k)**2)*(pr+perc) + CALL GR_TRANSFER(5._sp, ac_prcp(k), prr, ac_ct(k), ac_ht(k), & +& qr) + IF (0._sp .LT. prd + l) THEN + qd = prd + l + ELSE + qd = 0._sp + END IF + ac_qt(k) = qr + qd ! Transform from mm/dt to m3/s ac_qt(k) = ac_qt(k)*1e-3_sp*mesh%dx(row, col)*mesh%dy(row, col& & )/setup%dt @@ -15941,25 +15530,24 @@ SUBROUTINE GR4_MLP_ODE_TIME_STEP_D(setup, mesh, input_data, options, & END DO END DO output_layer_d = 0.0_4 - input_layer_d = 0.0_4 ! Forward MLP without OPENMP DO col=1,mesh%ncol DO row=1,mesh%nrow IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - input_layer_d(1) = ac_hp_d(k) - input_layer(1) = ac_hp(k) - input_layer_d(2) = ac_ht_d(k) - input_layer(2) = ac_ht(k) - input_layer_d(3) = pn_d(k) - input_layer(3) = pn(k) - input_layer_d(4) = en_d(k) - input_layer(4) = en(k) - CALL FORWARD_MLP_D(weight_1, weight_1_d, bias_1, bias_1_d, & -& weight_2, weight_2_d, bias_2, bias_2_d, & -& input_layer, input_layer_d, output_layer(:, k), & -& output_layer_d(:, k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + input_layer_d(:) = (/ac_hp_d(k), ac_ht_d(k), pn_d(k), en_d(k& +& )/) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP_D(weight_1, weight_1_d, bias_1, bias_1_d, & +& weight_2, weight_2_d, bias_2, bias_2_d, & +& input_layer, input_layer_d, output_layer(:, k)& +& , output_layer_d(:, k)) + ELSE + output_layer_d(:, k) = 0.0_4 + output_layer(:, k) = 0._sp + END IF END IF END DO END DO @@ -16082,20 +15670,19 @@ SUBROUTINE GR4_MLP_ODE_TIME_STEP_B(setup, mesh, input_data, options, & DO row=1,mesh%nrow IF (mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0) THEN - CALL PUSHCONTROL1B(0) + CALL PUSHCONTROL2B(0) ELSE k = mesh%rowcol_to_ind_ac(row, col) - CALL PUSHREAL4(input_layer(1)) - input_layer(1) = ac_hp(k) - CALL PUSHREAL4(input_layer(2)) - input_layer(2) = ac_ht(k) - CALL PUSHREAL4(input_layer(3)) - input_layer(3) = pn(k) - CALL PUSHREAL4(input_layer(4)) - input_layer(4) = en(k) - CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & -& input_layer, output_layer(:, k)) - CALL PUSHCONTROL1B(1) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + CALL PUSHREAL4ARRAY(input_layer, 4) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & +& input_layer, output_layer(:, k)) + CALL PUSHCONTROL2B(2) + ELSE + output_layer(:, k) = 0._sp + CALL PUSHCONTROL2B(1) + END IF END IF END DO END DO @@ -16157,29 +15744,26 @@ SUBROUTINE GR4_MLP_ODE_TIME_STEP_B(setup, mesh, input_data, options, & END DO END DO !$OMP END PARALLEL - input_layer_b = 0.0_4 DO col=mesh%ncol,1,-1 DO row=mesh%nrow,1,-1 - CALL POPCONTROL1B(branch) + CALL POPCONTROL2B(branch) IF (branch .NE. 0) THEN - k = mesh%rowcol_to_ind_ac(row, col) - CALL FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & -& weight_2, weight_2_b, bias_2, bias_2_b, & -& input_layer, input_layer_b, output_layer(:, k), & -& output_layer_b(:, k)) - output_layer_b(:, k) = 0.0_4 - CALL POPREAL4(input_layer(4)) - en_b(k) = en_b(k) + input_layer_b(4) - input_layer_b(4) = 0.0_4 - CALL POPREAL4(input_layer(3)) - pn_b(k) = pn_b(k) + input_layer_b(3) - input_layer_b(3) = 0.0_4 - CALL POPREAL4(input_layer(2)) - ac_ht_b(k) = ac_ht_b(k) + input_layer_b(2) - input_layer_b(2) = 0.0_4 - CALL POPREAL4(input_layer(1)) - ac_hp_b(k) = ac_hp_b(k) + input_layer_b(1) - input_layer_b(1) = 0.0_4 + IF (branch .EQ. 1) THEN + k = mesh%rowcol_to_ind_ac(row, col) + output_layer_b(:, k) = 0.0_4 + ELSE + k = mesh%rowcol_to_ind_ac(row, col) + CALL FORWARD_MLP_B(weight_1, weight_1_b, bias_1, bias_1_b, & +& weight_2, weight_2_b, bias_2, bias_2_b, & +& input_layer, input_layer_b, output_layer(:, k)& +& , output_layer_b(:, k)) + output_layer_b(:, k) = 0.0_4 + CALL POPREAL4ARRAY(input_layer, 4) + ac_hp_b(k) = ac_hp_b(k) + input_layer_b(1) + ac_ht_b(k) = ac_ht_b(k) + input_layer_b(2) + pn_b(k) = pn_b(k) + input_layer_b(3) + en_b(k) = en_b(k) + input_layer_b(4) + END IF END IF END DO END DO @@ -16273,12 +15857,13 @@ SUBROUTINE GR4_MLP_ODE_TIME_STEP(setup, mesh, input_data, options, & IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%& & local_active_cell(row, col) .EQ. 0)) THEN k = mesh%rowcol_to_ind_ac(row, col) - input_layer(1) = ac_hp(k) - input_layer(2) = ac_ht(k) - input_layer(3) = pn(k) - input_layer(4) = en(k) - CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & -& input_layer, output_layer(:, k)) + IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + CALL FORWARD_MLP(weight_1, bias_1, weight_2, bias_2, & +& input_layer, output_layer(:, k)) + ELSE + output_layer(:, k) = 0._sp + END IF END IF END DO END DO @@ -16361,9 +15946,9 @@ SUBROUTINE GR5_TIME_STEP_D(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION_D(ac_prcp(k), ac_prcp_d(k), ac_pet(k), & & ac_ci(k), ac_ci_d(k), ac_hi(k), ac_hi_d(k)& & , pn, pn_d, en, en_d) - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k& -& ), beta, ac_hp(k), ac_hp_d(k), pr, pr_d, perc& -& , perc_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_cp(k), ac_cp_d(k), beta, ac_hp(k& +& ), ac_hp_d(k), pr, pr_d, perc, perc_d) CALL GR_THRESHOLD_EXCHANGE_D(ac_kexc(k), ac_kexc_d(k), & & ac_aexc(k), ac_aexc_d(k), ac_ht(k), & & ac_ht_d(k), l, l_d) @@ -16432,6 +16017,8 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, returns, & REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd REAL(sp) :: pn_b, en_b, pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 INTEGER :: branch INTEGER :: chunk_start INTEGER :: chunk_end @@ -16462,8 +16049,8 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)& & , l) CALL PUSHCONTROL1B(1) @@ -16545,9 +16132,12 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, returns, & & ac_aexc(k), ac_aexc_b(k), ac_ht(k), & & ac_ht_b(k), l, l_b) CALL POPREAL4(ac_hp(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k& -& ), beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc& -& , perc_b) + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k), & +& beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc, & +& perc_b) CALL POPREAL4(ac_hi(k)) CALL POPREAL4(pn) CALL POPREAL4(en) @@ -16601,8 +16191,8 @@ SUBROUTINE GR5_TIME_STEP(setup, mesh, input_data, options, returns, & IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)& & , l) ELSE @@ -16689,9 +16279,9 @@ SUBROUTINE GR6_TIME_STEP_D(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION_D(ac_prcp(k), ac_prcp_d(k), ac_pet(k), & & ac_ci(k), ac_ci_d(k), ac_hi(k), ac_hi_d(k)& & , pn, pn_d, en, en_d) - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k& -& ), beta, ac_hp(k), ac_hp_d(k), pr, pr_d, perc& -& , perc_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_cp(k), ac_cp_d(k), beta, ac_hp(k& +& ), ac_hp_d(k), pr, pr_d, perc, perc_d) CALL GR_THRESHOLD_EXCHANGE_D(ac_kexc(k), ac_kexc_d(k), & & ac_aexc(k), ac_aexc_d(k), ac_ht(k), & & ac_ht_d(k), l, l_d) @@ -16767,6 +16357,8 @@ SUBROUTINE GR6_TIME_STEP_B(setup, mesh, input_data, options, returns, & REAL(sp) :: pn_b, en_b, pr_b, perc_b, l_b, prr_b, pre_b, prd_b, qr_b& & , qd_b, qe_b INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 REAL(sp) :: temp_b INTEGER :: branch INTEGER :: chunk_start @@ -16798,8 +16390,8 @@ SUBROUTINE GR6_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)& & , l) CALL PUSHCONTROL1B(1) @@ -16899,9 +16491,12 @@ SUBROUTINE GR6_TIME_STEP_B(setup, mesh, input_data, options, returns, & & ac_aexc(k), ac_aexc_b(k), ac_ht(k), & & ac_ht_b(k), l, l_b) CALL POPREAL4(ac_hp(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k& -& ), beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc& -& , perc_b) + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k), & +& beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc, & +& perc_b) CALL POPREAL4(ac_hi(k)) CALL POPREAL4(pn) CALL POPREAL4(en) @@ -16956,8 +16551,8 @@ SUBROUTINE GR6_TIME_STEP(setup, mesh, input_data, options, returns, & IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(& & k), pn, en) - CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), beta, & +& ac_hp(k), pr, perc) CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)& & , l) ELSE @@ -17049,9 +16644,9 @@ SUBROUTINE GRD_TIME_STEP_D(setup, mesh, input_data, options, returns, & END IF en_d = -ei_d en = ac_pet(k) - ei - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k& -& ), 1000._sp, ac_hp(k), ac_hp_d(k), pr, pr_d, & -& perc, perc_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_cp(k), ac_cp_d(k), 1000._sp, & +& ac_hp(k), ac_hp_d(k), pr, pr_d, perc, perc_d) ELSE pr = 0._sp perc = 0._sp @@ -17103,6 +16698,8 @@ SUBROUTINE GRD_TIME_STEP_B(setup, mesh, input_data, options, returns, & REAL(sp) :: ei_b, pn_b, en_b, pr_b, perc_b, prr_b, qr_b INTRINSIC MIN INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 INTEGER :: branch INTEGER :: chunk_start INTEGER :: chunk_end @@ -17143,8 +16740,8 @@ SUBROUTINE GRD_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL PUSHREAL4(en) en = ac_pet(k) - ei CALL PUSHREAL4(ac_hp(k)) - CALL GR_PRODUCTION(pn, en, ac_cp(k), 1000._sp, ac_hp(k), pr& -& , perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), 1000._sp& +& , ac_hp(k), pr, perc) CALL PUSHCONTROL1B(0) ELSE CALL PUSHCONTROL1B(1) @@ -17201,8 +16798,11 @@ SUBROUTINE GRD_TIME_STEP_B(setup, mesh, input_data, options, returns, & CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN CALL POPREAL4(ac_hp(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k& -& ), 1000._sp, ac_hp(k), ac_hp_b(k), pr, pr_b, & + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k), & +& 1000._sp, ac_hp(k), ac_hp_b(k), pr, pr_b, & & perc, perc_b) CALL POPREAL4(en) ei_b = -en_b @@ -17272,8 +16872,8 @@ SUBROUTINE GRD_TIME_STEP(setup, mesh, input_data, options, returns, & pn = 0._sp END IF en = ac_pet(k) - ei - CALL GR_PRODUCTION(pn, en, ac_cp(k), 1000._sp, ac_hp(k), pr& -& , perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_cp(k), 1000._sp& +& , ac_hp(k), pr, perc) ELSE pr = 0._sp perc = 0._sp @@ -17358,9 +16958,9 @@ SUBROUTINE LOIEAU_TIME_STEP_D(setup, mesh, input_data, options, & END IF en_d = -ei_d en = ac_pet(k) - ei - CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_ca(k), ac_ca_d(k& -& ), beta, ac_ha(k), ac_ha_d(k), pr, pr_d, perc& -& , perc_d) + CALL GR_PRODUCTION_D(0._sp, 0.0_4, 0._sp, 0.0_4, pn, pn_d, & +& en, en_d, ac_ca(k), ac_ca_d(k), beta, ac_ha(k& +& ), ac_ha_d(k), pr, pr_d, perc, perc_d) ELSE pr = 0._sp perc = 0._sp @@ -17421,6 +17021,8 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, & REAL(sp) :: ei_b, pn_b, en_b, pr_b, perc_b, prr_b, prd_b, qr_b, qd_b INTRINSIC MIN INTRINSIC MAX + REAL(sp) :: dummydiff_b + REAL(sp) :: dummydiff_b0 INTEGER :: branch INTEGER :: chunk_start INTEGER :: chunk_end @@ -17463,8 +17065,8 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, & CALL PUSHREAL4(en) en = ac_pet(k) - ei CALL PUSHREAL4(ac_ha(k)) - CALL GR_PRODUCTION(pn, en, ac_ca(k), beta, ac_ha(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_ca(k), beta, & +& ac_ha(k), pr, perc) CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) @@ -17550,9 +17152,12 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, & CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN CALL POPREAL4(ac_ha(k)) - CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_ca(k), ac_ca_b(k& -& ), beta, ac_ha(k), ac_ha_b(k), pr, pr_b, perc& -& , perc_b) + pn_b = 0.0_4 + en_b = 0.0_4 + CALL GR_PRODUCTION_B(0._sp, dummydiff_b, 0._sp, dummydiff_b0& +& , pn, pn_b, en, en_b, ac_ca(k), ac_ca_b(k), & +& beta, ac_ha(k), ac_ha_b(k), pr, pr_b, perc, & +& perc_b) CALL POPREAL4(en) ei_b = -en_b CALL POPCONTROL1B(branch) @@ -17623,8 +17228,8 @@ SUBROUTINE LOIEAU_TIME_STEP(setup, mesh, input_data, options, returns& pn = 0._sp END IF en = ac_pet(k) - ei - CALL GR_PRODUCTION(pn, en, ac_ca(k), beta, ac_ha(k), pr, & -& perc) + CALL GR_PRODUCTION(0._sp, 0._sp, pn, en, ac_ca(k), beta, & +& ac_ha(k), pr, perc) ELSE pr = 0._sp perc = 0._sp diff --git a/smash/fcore/operator/md_gr_operator.f90 b/smash/fcore/operator/md_gr_operator.f90 index f141a484..c418d013 100644 --- a/smash/fcore/operator/md_gr_operator.f90 +++ b/smash/fcore/operator/md_gr_operator.f90 @@ -8,7 +8,6 @@ !% - gr_exchange !% - gr_threshold_exchange !% - gr_transfer -!% - gr_production_transfer_mlp_alg !% - gr_production_transfer_ode !% - gr_production_transfer_mlp_ode !% - gr4_time_step @@ -56,11 +55,11 @@ subroutine gr_interception(prcp, pet, ci, hi, pn, en) end subroutine gr_interception - subroutine gr_production(pn, en, cp, beta, hp, pr, perc) + subroutine gr_production(fq_ps, fq_es, pn, en, cp, beta, hp, pr, perc) implicit none - real(sp), intent(in) :: pn, en, cp, beta + real(sp), intent(in) :: fq_ps, fq_es, pn, en, cp, beta real(sp), intent(inout) :: hp real(sp), intent(out) :: pr, perc @@ -71,9 +70,11 @@ subroutine gr_production(pn, en, cp, beta, hp, pr, perc) ps = cp*(1._sp - hp*hp)*tanh(pn*inv_cp)/ & & (1._sp + hp*tanh(pn*inv_cp)) + ps = (1._sp + fq_ps)*ps es = (hp*cp)*(2._sp - hp)*tanh(en*inv_cp)/ & & (1._sp + (1._sp - hp)*tanh(en*inv_cp)) + es = (1._sp + fq_es)*es hp_imd = hp + (ps - es)*inv_cp @@ -89,15 +90,15 @@ subroutine gr_production(pn, en, cp, beta, hp, pr, perc) end subroutine gr_production - subroutine gr_exchange(kexc, ht, l) + subroutine gr_exchange(fq_l, kexc, ht, l) implicit none - real(sp), intent(in) :: kexc + real(sp), intent(in) :: fq_l, kexc real(sp), intent(inout) :: ht real(sp), intent(out) :: l - l = kexc*ht**3.5_sp + l = (1._sp + fq_l)*kexc*ht**3.5_sp end subroutine gr_exchange @@ -167,74 +168,6 @@ subroutine gr_exponential_transfer(pre, be, he, qe) end subroutine gr_exponential_transfer - subroutine gr_production_transfer_mlp_alg(f_q, pn, en, cp, beta, ct, kexc, n, prcp, hp, ht, q, & - & pr, perc, l, prr, prd, qr, qd) - !% Integrate MLP in stepwise aprroximation method - - implicit none - - real(sp), dimension(4), intent(in) :: f_q ! fixed NN output size - real(sp), intent(in) :: pn, en, cp, beta, ct, kexc, n, prcp - real(sp), intent(inout) :: hp, ht, q - real(sp), intent(out) :: pr, perc, l, prr, prd, qr, qd - - real(sp) :: inv_cp, ps, es, hp_imd, pr_imd, ht_imd, nm1, d1pnm1 - - !% Production - - inv_cp = 1._sp/cp - pr = 0._sp - - ps = cp*(1._sp - hp*hp)*tanh(pn*inv_cp)/(1._sp + hp*tanh(pn*inv_cp)) - ps = (1._sp + f_q(1))*ps - - es = (hp*cp)*(2._sp - hp)*tanh(en*inv_cp)/(1._sp + (1._sp - hp)*tanh(en*inv_cp)) - es = (1._sp + f_q(2))*es - - hp_imd = hp + (ps - es)*inv_cp - - if (pn .gt. 0) then - - pr = pn - (hp_imd - hp)*cp - - end if - - perc = (hp_imd*cp)*(1._sp - (1._sp + (hp_imd/beta)**4)**(-0.25_sp)) - - hp = hp_imd - perc*inv_cp - - !% Transfer - - nm1 = n - 1._sp - d1pnm1 = 1._sp/nm1 - - l = (1._sp + f_q(4))*kexc*ht**3.5_sp - - prr = (0.9_sp*(1._sp - f_q(3)**2))*(pr + perc) + l - prd = (0.1_sp + 0.9_sp*f_q(3)**2)*(pr + perc) - - if (prcp .lt. 0._sp) then - - pr_imd = ((ht*ct)**(-nm1) - ct**(-nm1))**(-d1pnm1) - (ht*ct) - - else - - pr_imd = prr - - end if - - ht_imd = max(1.e-6_sp, ht + pr_imd/ct) - - ht = (((ht_imd*ct)**(-nm1) + ct**(-nm1))**(-d1pnm1))/ct - - qd = max(0._sp, prd + l) - - qr = (ht_imd - ht)*ct - - q = qr + qd - - end subroutine gr_production_transfer_mlp_alg - subroutine gr_production_transfer_ode(pn, en, cp, ct, kexc, hp, ht, q, l) !% Solve state-space ODE system with implicit Euler @@ -305,12 +238,12 @@ subroutine gr_production_transfer_ode(pn, en, cp, ct, kexc, hp, ht, q, l) end subroutine gr_production_transfer_ode - subroutine gr_production_transfer_mlp_ode(f_q, pn, en, cp, ct, kexc, hp, ht, q, l) + subroutine gr_production_transfer_mlp_ode(fq, pn, en, cp, ct, kexc, hp, ht, q, l) !% Solve state-space ODE system with explicit Euler and MLP implicit none - real(sp), dimension(5), intent(in) :: f_q ! fixed NN output size + real(sp), dimension(5), intent(in) :: fq ! fixed NN output size real(sp), intent(in) :: pn, en, cp, ct, kexc real(sp), intent(inout) :: hp, ht, q real(sp), intent(out) :: l @@ -326,17 +259,17 @@ subroutine gr_production_transfer_mlp_ode(f_q, pn, en, cp, ct, kexc, hp, ht, q, !do i = 1, n_subtimesteps - fhp = (1._sp + f_q(1))*pn*(1._sp - hp**2) & - & - (1._sp + f_q(2))*en*hp*(2._sp - hp) + fhp = (1._sp + fq(1))*pn*(1._sp - hp**2) & + & - (1._sp + fq(2))*en*hp*(2._sp - hp) hp = hp + dt*fhp*inv_cp if (hp .le. 0._sp) hp = 1.e-6_sp if (hp .ge. 1._sp) hp = 1._sp - 1.e-6_sp - fht = (0.9_sp*(1._sp - f_q(3)**2))*(1._sp + f_q(1))*pn*hp**2 & - & + (1._sp + f_q(4))*kexc*ht**3.5_sp & - & - (1._sp + f_q(5))*ct*ht**5 + fht = (0.9_sp*(1._sp - fq(3)**2))*(1._sp + fq(1))*pn*hp**2 & + & + (1._sp + fq(4))*kexc*ht**3.5_sp & + & - (1._sp + fq(5))*ct*ht**5 ht = ht + dt*fht/ct @@ -345,10 +278,10 @@ subroutine gr_production_transfer_mlp_ode(f_q, pn, en, cp, ct, kexc, hp, ht, q, !end do - l = (1._sp + f_q(4))*kexc*ht**3.5_sp + l = (1._sp + fq(4))*kexc*ht**3.5_sp - q = (1._sp + f_q(5))*ct*ht**5 & - & + (0.1_sp + 0.9_sp*f_q(3)**2)*(1._sp + f_q(1))*pn*hp**2 + l + q = (1._sp + fq(5))*ct*ht**5 & + & + (0.1_sp + 0.9_sp*fq(3)**2)*(1._sp + fq(1))*pn*hp**2 + l end subroutine gr_production_transfer_mlp_ode @@ -397,9 +330,9 @@ subroutine gr4_time_step(setup, mesh, input_data, options, returns, time_step, a call gr_interception(ac_prcp(k), ac_pet(k), ac_ci(k), & & ac_hi(k), pn, en) - call gr_production(pn, en, ac_cp(k), beta, ac_hp(k), pr, perc) + call gr_production(0._sp, 0._sp, pn, en, ac_cp(k), beta, ac_hp(k), pr, perc) - call gr_exchange(ac_kexc(k), ac_ht(k), l) + call gr_exchange(0._sp, ac_kexc(k), ac_ht(k), l) else @@ -519,12 +452,15 @@ subroutine gr4_mlp_alg_time_step(setup, mesh, input_data, options, returns, time k = mesh%rowcol_to_ind_ac(row, col) - input_layer(1) = ac_hp(k) - input_layer(2) = ac_ht(k) - input_layer(3) = pn(k) - input_layer(4) = en(k) + if (ac_prcp(k) .ge. 0._sp .and. ac_pet(k) .ge. 0._sp) then + + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + call forward_mlp(weight_1, bias_1, weight_2, bias_2, input_layer, output_layer(:, k)) - call forward_mlp(weight_1, bias_1, weight_2, bias_2, input_layer, output_layer(:, k)) + else + output_layer(:, k) = 0._sp + + end if end do end do @@ -532,7 +468,7 @@ subroutine gr4_mlp_alg_time_step(setup, mesh, input_data, options, returns, time ! Production and transfer with OPENMP #ifdef _OPENMP !$OMP parallel do schedule(static) num_threads(options%comm%ncpu) & - !$OMP& shared(setup, mesh, returns, output_layer, ac_prcp, & + !$OMP& shared(setup, mesh, returns, output_layer, ac_prcp, ac_pet, & !$OMP& ac_cp, beta, ac_ct, ac_kexc, ac_hp, ac_ht, ac_qt, pn, en) & !$OMP& private(row, col, k, time_step_returns, pr, perc, l, prr, prd, qr, qd) #endif @@ -543,9 +479,29 @@ subroutine gr4_mlp_alg_time_step(setup, mesh, input_data, options, returns, time k = mesh%rowcol_to_ind_ac(row, col) - call gr_production_transfer_mlp_alg(output_layer(:, k), pn(k), en(k), & - & ac_cp(k), beta, ac_ct(k), ac_kexc(k), 5._sp, ac_prcp(k), ac_hp(k), & - & ac_ht(k), ac_qt(k), pr, perc, l, prr, prd, qr, qd) + if (ac_prcp(k) .ge. 0._sp .and. ac_pet(k) .ge. 0._sp) then + + call gr_production(output_layer(1, k), output_layer(2, k), pn(k), en(k), ac_cp(k), & + & beta, ac_hp(k), pr, perc) + + call gr_exchange(output_layer(4, k), ac_kexc(k), ac_ht(k), l) + + else + + pr = 0._sp + perc = 0._sp + l = 0._sp + + end if + + prr = (0.9_sp*(1._sp - output_layer(3, k)**2))*(pr + perc) + l + prd = (0.1_sp + 0.9_sp*output_layer(3, k)**2)*(pr + perc) + + call gr_transfer(5._sp, ac_prcp(k), prr, ac_ct(k), ac_ht(k), qr) + + qd = max(0._sp, prd + l) + + ac_qt(k) = qr + qd ! Transform from mm/dt to m3/s ac_qt(k) = ac_qt(k)*1e-3_sp*mesh%dx(row, col)*mesh%dy(row, col)/setup%dt @@ -740,12 +696,15 @@ subroutine gr4_mlp_ode_time_step(setup, mesh, input_data, options, returns, time k = mesh%rowcol_to_ind_ac(row, col) - input_layer(1) = ac_hp(k) - input_layer(2) = ac_ht(k) - input_layer(3) = pn(k) - input_layer(4) = en(k) + if (ac_prcp(k) .ge. 0._sp .and. ac_pet(k) .ge. 0._sp) then - call forward_mlp(weight_1, bias_1, weight_2, bias_2, input_layer, output_layer(:, k)) + input_layer(:) = (/ac_hp(k), ac_ht(k), pn(k), en(k)/) + call forward_mlp(weight_1, bias_1, weight_2, bias_2, input_layer, output_layer(:, k)) + + else + output_layer(:, k) = 0._sp + + end if end do end do @@ -842,7 +801,7 @@ subroutine gr5_time_step(setup, mesh, input_data, options, returns, time_step, a call gr_interception(ac_prcp(k), ac_pet(k), ac_ci(k), & & ac_hi(k), pn, en) - call gr_production(pn, en, ac_cp(k), beta, ac_hp(k), pr, perc) + call gr_production(0._sp, 0._sp, pn, en, ac_cp(k), beta, ac_hp(k), pr, perc) call gr_threshold_exchange(ac_kexc(k), ac_aexc(k), ac_ht(k), l) @@ -937,7 +896,7 @@ subroutine gr6_time_step(setup, mesh, input_data, options, returns, time_step, a call gr_interception(ac_prcp(k), ac_pet(k), ac_ci(k), & & ac_hi(k), pn, en) - call gr_production(pn, en, ac_cp(k), beta, ac_hp(k), pr, perc) + call gr_production(0._sp, 0._sp, pn, en, ac_cp(k), beta, ac_hp(k), pr, perc) call gr_threshold_exchange(ac_kexc(k), ac_aexc(k), ac_ht(k), l) @@ -1034,7 +993,7 @@ subroutine grd_time_step(setup, mesh, input_data, options, returns, time_step, a en = ac_pet(k) - ei - call gr_production(pn, en, ac_cp(k), 1000._sp, ac_hp(k), pr, perc) + call gr_production(0._sp, 0._sp, pn, en, ac_cp(k), 1000._sp, ac_hp(k), pr, perc) else @@ -1125,7 +1084,7 @@ subroutine loieau_time_step(setup, mesh, input_data, options, returns, time_step en = ac_pet(k) - ei - call gr_production(pn, en, ac_ca(k), beta, ac_ha(k), pr, perc) + call gr_production(0._sp, 0._sp, pn, en, ac_ca(k), beta, ac_ha(k), pr, perc) else