1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculates the exchange energy based on the Becke-Roussel exchange 8!> hole. Takes advantage of an analytical representation of the hole 9!> in order to avoid solving a non-linear equation by means of Newton- 10!> Raphson algorithm 11!> \note 12!> \par History 13!> 12.2008 created [mguidon] 14!> \author mguidon 15! ************************************************************************************************** 16MODULE xc_xbecke_roussel 17 USE bibliography, ONLY: BeckeRoussel1989,& 18 Proynov2007,& 19 cite_reference 20 USE input_section_types, ONLY: section_vals_type,& 21 section_vals_val_get 22 USE kinds, ONLY: dp 23 USE mathconstants, ONLY: pi 24 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 25 xc_dset_get_derivative 26 USE xc_derivative_types, ONLY: xc_derivative_get,& 27 xc_derivative_type 28 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 29 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 30 xc_rho_set_type 31#include "../base/base_uses.f90" 32 33 IMPLICIT NONE 34 PRIVATE 35 36 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbecke_roussel' 38 39 REAL(dp), PARAMETER, PRIVATE :: br_a1 = 1.5255251812009530_dp, & 40 br_a2 = 0.4576575543602858_dp, & 41 br_a3 = 0.4292036732051034_dp, & 42 br_c0 = 0.7566445420735584_dp, & 43 br_c1 = -2.6363977871370960_dp, & 44 br_c2 = 5.4745159964232880_dp, & 45 br_c3 = -12.657308127108290_dp, & 46 br_c4 = 4.1250584725121360_dp, & 47 br_c5 = -30.425133957163840_dp, & 48 br_b0 = 0.4771976183772063_dp, & 49 br_b1 = -1.7799813494556270_dp, & 50 br_b2 = 3.8433841862302150_dp, & 51 br_b3 = -9.5912050880518490_dp, & 52 br_b4 = 2.1730180285916720_dp, & 53 br_b5 = -30.425133851603660_dp, & 54 br_d0 = 0.00004435009886795587_dp, & 55 br_d1 = 0.58128653604457910_dp, & 56 br_d2 = 66.742764515940610_dp, & 57 br_d3 = 434.26780897229770_dp, & 58 br_d4 = 824.7765766052239000_dp, & 59 br_d5 = 1657.9652731582120_dp, & 60 br_e0 = 0.00003347285060926091_dp, & 61 br_e1 = 0.47917931023971350_dp, & 62 br_e2 = 62.392268338574240_dp, & 63 br_e3 = 463.14816427938120_dp, & 64 br_e4 = 785.2360350104029000_dp, & 65 br_e5 = 1657.962968223273000000_dp, & 66 br_BB = 2.085749716493756_dp 67 68 PUBLIC :: xbecke_roussel_lda_info, xbecke_roussel_lsd_info, xbecke_roussel_lda_eval, & 69 xbecke_roussel_lsd_eval, x_br_lsd_y_lte_0, x_br_lsd_y_gt_0, x_br_lsd_y_lte_0_cutoff, & 70 x_br_lsd_y_gt_0_cutoff 71CONTAINS 72 73! ************************************************************************************************** 74!> \brief return various information on the functional 75!> \param reference string with the reference of the actual functional 76!> \param shortform string with the shortform of the functional name 77!> \param needs the components needed by this functional are set to 78!> true (does not set the unneeded components to false) 79!> \param max_deriv controls the number of derivatives 80!> \par History 81!> 12.2008 created [mguidon] 82!> \author mguidon 83! ************************************************************************************************** 84 SUBROUTINE xbecke_roussel_lda_info(reference, shortform, needs, max_deriv) 85 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 86 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs 87 INTEGER, INTENT(out), OPTIONAL :: max_deriv 88 89 CHARACTER(len=*), PARAMETER :: routineN = 'xbecke_roussel_lda_info', & 90 routineP = moduleN//':'//routineN 91 92 CALL cite_reference(BeckeRoussel1989) 93 CALL cite_reference(Proynov2007) 94 IF (PRESENT(reference)) THEN 95 reference = "A.D. Becke, M.R. Roussel, "// & 96 "Phys. Rev. A, vol. 39, n 8, pp. 3761-3767, (1989)"// & 97 "{spin unpolarized}" 98 END IF 99 IF (PRESENT(shortform)) THEN 100 shortform = "Becke-Roussel exchange hole (spin unpolarized)" 101 END IF 102 IF (PRESENT(needs)) THEN 103 needs%rho = .TRUE. 104 needs%norm_drho = .TRUE. 105 needs%tau = .TRUE. 106 needs%laplace_rho = .TRUE. 107 END IF 108 109 IF (PRESENT(max_deriv)) max_deriv = 1 110 111 END SUBROUTINE xbecke_roussel_lda_info 112 113! ************************************************************************************************** 114!> \brief return various information on the functional 115!> \param reference string with the reference of the actual functional 116!> \param shortform string with the shortform of the functional name 117!> \param needs the components needed by this functional are set to 118!> true (does not set the unneeded components to false) 119!> \param max_deriv ... 120!> \par History 121!> 12.2008 created [mguidon] 122!> \author mguidon 123! ************************************************************************************************** 124 SUBROUTINE xbecke_roussel_lsd_info(reference, shortform, needs, max_deriv) 125 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 126 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs 127 INTEGER, INTENT(out), OPTIONAL :: max_deriv 128 129 CHARACTER(len=*), PARAMETER :: routineN = 'xbecke_roussel_lsd_info', & 130 routineP = moduleN//':'//routineN 131 132 CALL cite_reference(BeckeRoussel1989) 133 CALL cite_reference(Proynov2007) 134 IF (PRESENT(reference)) THEN 135 reference = "A.D. Becke, M.R. Roussel, "// & 136 "Phys. Rev. A, vol. 39, n 8, pp. 3761-3767, (1989)"// & 137 "{spin polarized}" 138 END IF 139 IF (PRESENT(shortform)) THEN 140 shortform = "Becke-Roussel exchange hole (spin polarized)" 141 END IF 142 IF (PRESENT(needs)) THEN 143 needs%rho_spin = .TRUE. 144 needs%norm_drho_spin = .TRUE. 145 needs%tau_spin = .TRUE. 146 needs%laplace_rho_spin = .TRUE. 147 END IF 148 IF (PRESENT(max_deriv)) max_deriv = 1 149 150 END SUBROUTINE xbecke_roussel_lsd_info 151 152! ************************************************************************************************** 153!> \brief evaluates the Becke Roussel exchange functional for lda 154!> \param rho_set the density where you want to evaluate the functional 155!> \param deriv_set place where to store the functional derivatives (they are 156!> added to the derivatives) 157!> \param grad_deriv degree of the derivative that should be evaluated, 158!> if positive all the derivatives up to the given degree are evaluated, 159!> if negative only the given degree is calculated 160!> \param br_params parameters for the becke roussel functional 161!> \par History 162!> 12.2008 created [mguidon] 163!> \author mguidon 164! ************************************************************************************************** 165 SUBROUTINE xbecke_roussel_lda_eval(rho_set, deriv_set, grad_deriv, br_params) 166 TYPE(xc_rho_set_type), POINTER :: rho_set 167 TYPE(xc_derivative_set_type), POINTER :: deriv_set 168 INTEGER, INTENT(in) :: grad_deriv 169 TYPE(section_vals_type), POINTER :: br_params 170 171 CHARACTER(len=*), PARAMETER :: routineN = 'xbecke_roussel_lda_eval', & 172 routineP = moduleN//':'//routineN 173 174 INTEGER :: handle, npoints 175 INTEGER, DIMENSION(:, :), POINTER :: bo 176 REAL(dp) :: gamma, R, sx 177 REAL(kind=dp) :: epsilon_rho 178 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rho, e_ndrho, & 179 e_rho, e_tau, laplace_rho, norm_drho, & 180 rho, tau 181 TYPE(xc_derivative_type), POINTER :: deriv 182 183 CALL timeset(routineN, handle) 184 185 NULLIFY (bo) 186 187 CPASSERT(ASSOCIATED(rho_set)) 188 CPASSERT(rho_set%ref_count > 0) 189 CPASSERT(ASSOCIATED(deriv_set)) 190 CPASSERT(deriv_set%ref_count > 0) 191 CALL xc_rho_set_get(rho_set, rho=rho, norm_drho=norm_drho, & 192 tau=tau, laplace_rho=laplace_rho, local_bounds=bo, & 193 rho_cutoff=epsilon_rho) 194 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 195 196 dummy => rho 197 198 e_0 => dummy 199 e_rho => dummy 200 e_ndrho => dummy 201 e_tau => dummy 202 e_laplace_rho => dummy 203 204 IF (grad_deriv >= 0) THEN 205 deriv => xc_dset_get_derivative(deriv_set, "", & 206 allocate_deriv=.TRUE.) 207 CALL xc_derivative_get(deriv, deriv_data=e_0) 208 END IF 209 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 210 deriv => xc_dset_get_derivative(deriv_set, "(rho)", & 211 allocate_deriv=.TRUE.) 212 CALL xc_derivative_get(deriv, deriv_data=e_rho) 213 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 214 allocate_deriv=.TRUE.) 215 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 216 deriv => xc_dset_get_derivative(deriv_set, "(tau)", & 217 allocate_deriv=.TRUE.) 218 CALL xc_derivative_get(deriv, deriv_data=e_tau) 219 deriv => xc_dset_get_derivative(deriv_set, "(laplace_rho)", & 220 allocate_deriv=.TRUE.) 221 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho) 222 END IF 223 IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN 224 CPABORT("derivatives bigger than 1 not implemented") 225 END IF 226 227 CALL section_vals_val_get(br_params, "scale_x", r_val=sx) 228 CALL section_vals_val_get(br_params, "CUTOFF_RADIUS", r_val=R) 229 CALL section_vals_val_get(br_params, "GAMMA", r_val=gamma) 230 231!$OMP PARALLEL DEFAULT(NONE) & 232!$OMP SHARED(rho, norm_drho, laplace_rho, tau, e_0, e_rho) & 233!$OMP SHARED(e_ndrho, e_tau, e_laplace_rho, grad_deriv) & 234!$OMP SHARED(npoints, epsilon_rho) & 235!$OMP SHARED(sx, r, gamma) 236 237 CALL xbecke_roussel_lda_calc(rho=rho, norm_drho=norm_drho, & 238 laplace_rho=laplace_rho, tau=tau, e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, & 239 e_tau=e_tau, e_laplace_rho=e_laplace_rho, grad_deriv=grad_deriv, & 240 npoints=npoints, epsilon_rho=epsilon_rho, & 241 sx=sx, R=R, gamma=gamma) 242 243!$OMP END PARALLEL 244 245 CALL timestop(handle) 246 END SUBROUTINE xbecke_roussel_lda_eval 247 248! ************************************************************************************************** 249!> \brief Precalculates which branch of the code has to be taken 250!> There are two main branches, one for a truncated potential and a cutoff 251!> radius, the other for the full coulomb interaction. In the end, the code 252!> for the lsd part will be called and the lda part is obtained via spin 253!> scaling relations 254!> \param rho grid values 255!> \param norm_drho grid values 256!> \param laplace_rho grid values 257!> \param tau grid values 258!> \param e_0 energies and derivatives 259!> \param e_rho energies and derivatives 260!> \param e_ndrho energies and derivatives 261!> \param e_tau energies and derivatives 262!> \param e_laplace_rho energies and derivatives 263!> \param grad_deriv degree of the derivative that should be evaluated, 264!> if positive all the derivatives up to the given degree are evaluated, 265!> if negative only the given degree is calculated 266!> \param npoints size of the grids 267!> \param epsilon_rho cutoffs 268!> \param sx scales the exchange potential and energies 269!> \param R cutoff Radius for truncated case 270!> \param gamma parameter from original publication, usually set to 1 271!> \par History 272!> 12.2008 created [mguidon] 273!> \author mguidon 274! ************************************************************************************************** 275 SUBROUTINE xbecke_roussel_lda_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, & 276 e_ndrho, e_tau, e_laplace_rho, grad_deriv, npoints, & 277 epsilon_rho, sx, R, gamma) 278 279 INTEGER, INTENT(in) :: npoints, grad_deriv 280 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_laplace_rho, e_tau, e_ndrho, e_rho, e_0 281 REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: tau, laplace_rho, norm_drho, rho 282 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R, gamma 283 284 CHARACTER(len=*), PARAMETER :: routineN = 'xbecke_roussel_lda_calc', & 285 routineP = moduleN//':'//routineN 286 287 INTEGER :: ip 288 REAL(dp) :: e_diff, e_old, my_laplace_rho, my_ndrho, & 289 my_rho, my_tau, t1, t15, t16, t2, t3, & 290 t4, t5, t8, t9, yval 291 292! Precalculate y, in order to chose the correct branch afterwards 293 294!$OMP DO 295 296 DO ip = 1, npoints 297 my_rho = 0.5_dp*MAX(rho(ip), 0.0_dp) 298 IF (my_rho > epsilon_rho) THEN 299 my_ndrho = 0.5_dp*MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp) 300 my_tau = 0.5_dp*MAX(EPSILON(0.0_dp)*1.e4_dp, tau(ip)) 301 my_laplace_rho = 0.5_dp*laplace_rho(ip) 302 303 t1 = pi**(0.1e1_dp/0.3e1_dp) 304 t2 = t1**2 305 t3 = my_rho**(0.1e1_dp/0.3e1_dp) 306 t4 = t3**2 307 t5 = t4*my_rho 308 t8 = my_ndrho**2 309 t9 = 0.1e1_dp/my_rho 310 ! *** CP2K defines tau in a different way as compared to Becke !!! 311 t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp 312 t16 = 0.1e1_dp/t15 313 yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16 314 IF (R == 0.0_dp) THEN 315 IF (yval <= 0.0_dp) THEN 316 e_old = e_0(ip) 317 CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), & 318 e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), & 319 sx, gamma, grad_deriv) 320 ! VERY UGLY HACK e_0 has to multiplied by the factor 2 321 e_diff = e_0(ip) - e_old 322 e_0(ip) = e_0(ip) + e_diff 323 ELSE 324 e_old = e_0(ip) 325 CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), & 326 e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), & 327 sx, gamma, grad_deriv) 328 ! VERY UGLY HACK e_0 has to multiplied by the factor 2 329 e_diff = e_0(ip) - e_old 330 e_0(ip) = e_0(ip) + e_diff 331 END IF 332 ELSE 333 IF (yval <= 0.0_dp) THEN 334 e_old = e_0(ip) 335 CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), & 336 e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), & 337 sx, R, gamma, grad_deriv) 338 ! VERY UGLY HACK e_0 has to multiplied by the factor 2 339 e_diff = e_0(ip) - e_old 340 e_0(ip) = e_0(ip) + e_diff 341 ELSE 342 e_old = e_0(ip) 343 CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), & 344 e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), & 345 sx, R, gamma, grad_deriv) 346 ! VERY UGLY HACK e_0 has to multiplied by the factor 2 347 e_diff = e_0(ip) - e_old 348 e_0(ip) = e_0(ip) + e_diff 349 END IF 350 END IF 351 END IF 352 END DO 353 354!$OMP END DO 355 356 END SUBROUTINE xbecke_roussel_lda_calc 357 358! ************************************************************************************************** 359!> \brief evaluates the Becke Roussel exchange functional for lda 360!> \param rho_set the density where you want to evaluate the functional 361!> \param deriv_set place where to store the functional derivatives (they are 362!> added to the derivatives) 363!> \param grad_deriv degree of the derivative that should be evaluated, 364!> if positive all the derivatives up to the given degree are evaluated, 365!> if negative only the given degree is calculated 366!> \param br_params parameters for the becke roussel functional 367!> \par History 368!> 12.2008 created [mguidon] 369!> \author mguidon 370! ************************************************************************************************** 371 SUBROUTINE xbecke_roussel_lsd_eval(rho_set, deriv_set, grad_deriv, br_params) 372 TYPE(xc_rho_set_type), POINTER :: rho_set 373 TYPE(xc_derivative_set_type), POINTER :: deriv_set 374 INTEGER, INTENT(in) :: grad_deriv 375 TYPE(section_vals_type), POINTER :: br_params 376 377 CHARACTER(len=*), PARAMETER :: routineN = 'xbecke_roussel_lsd_eval', & 378 routineP = moduleN//':'//routineN 379 380 INTEGER :: handle, npoints 381 INTEGER, DIMENSION(:, :), POINTER :: bo 382 REAL(dp) :: gamma, R, sx 383 REAL(kind=dp) :: epsilon_rho 384 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, e_laplace_rhob, & 385 e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, e_tau_a, e_tau_b, laplace_rhoa, laplace_rhob, & 386 norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b 387 TYPE(xc_derivative_type), POINTER :: deriv 388 389 CALL timeset(routineN, handle) 390 391 NULLIFY (bo) 392 393 CPASSERT(ASSOCIATED(rho_set)) 394 CPASSERT(rho_set%ref_count > 0) 395 CPASSERT(ASSOCIATED(deriv_set)) 396 CPASSERT(deriv_set%ref_count > 0) 397 CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, & 398 norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, laplace_rhoa=laplace_rhoa, & 399 laplace_rhob=laplace_rhob, local_bounds=bo, & 400 rho_cutoff=epsilon_rho) 401 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 402 403 dummy => rhoa 404 405 e_0 => dummy 406 e_rhoa => dummy 407 e_rhob => dummy 408 e_ndrhoa => dummy 409 e_ndrhob => dummy 410 e_tau_a => dummy 411 e_tau_b => dummy 412 e_laplace_rhoa => dummy 413 e_laplace_rhob => dummy 414 415 IF (grad_deriv >= 0) THEN 416 deriv => xc_dset_get_derivative(deriv_set, "", & 417 allocate_deriv=.TRUE.) 418 CALL xc_derivative_get(deriv, deriv_data=e_0) 419 END IF 420 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 421 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)", & 422 allocate_deriv=.TRUE.) 423 CALL xc_derivative_get(deriv, deriv_data=e_rhoa) 424 deriv => xc_dset_get_derivative(deriv_set, "(rhob)", & 425 allocate_deriv=.TRUE.) 426 CALL xc_derivative_get(deriv, deriv_data=e_rhob) 427 428 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)", & 429 allocate_deriv=.TRUE.) 430 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa) 431 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)", & 432 allocate_deriv=.TRUE.) 433 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob) 434 435 deriv => xc_dset_get_derivative(deriv_set, "(tau_a)", & 436 allocate_deriv=.TRUE.) 437 CALL xc_derivative_get(deriv, deriv_data=e_tau_a) 438 deriv => xc_dset_get_derivative(deriv_set, "(tau_b)", & 439 allocate_deriv=.TRUE.) 440 CALL xc_derivative_get(deriv, deriv_data=e_tau_b) 441 442 deriv => xc_dset_get_derivative(deriv_set, "(laplace_rhoa)", & 443 allocate_deriv=.TRUE.) 444 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa) 445 deriv => xc_dset_get_derivative(deriv_set, "(laplace_rhob)", & 446 allocate_deriv=.TRUE.) 447 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob) 448 END IF 449 IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN 450 CPABORT("derivatives bigger than 1 not implemented") 451 END IF 452 453 CALL section_vals_val_get(br_params, "scale_x", r_val=sx) 454 CALL section_vals_val_get(br_params, "CUTOFF_RADIUS", r_val=R) 455 CALL section_vals_val_get(br_params, "GAMMA", r_val=gamma) 456 457!$OMP PARALLEL DEFAULT (NONE) & 458!$OMP SHARED(rhoa, norm_drhoa, laplace_rhoa, tau_a, e_0) & 459!$OMP SHARED(e_rhoa, e_ndrhoa, e_tau_a, e_laplace_rhoa) & 460!$OMP SHARED(grad_deriv, npoints, epsilon_rho) & 461!$OMP SHARED(sx, r, gamma) & 462!$OMP SHARED(rhob, norm_drhob, laplace_rhob, tau_b, e_rhob) & 463!$OMP SHARED(e_ndrhob, e_tau_b, e_laplace_rhob) 464 465 CALL xbecke_roussel_lsd_calc(rho=rhoa, norm_drho=norm_drhoa, & 466 laplace_rho=laplace_rhoa, tau=tau_a, e_0=e_0, e_rho=e_rhoa, e_ndrho=e_ndrhoa, & 467 e_tau=e_tau_a, e_laplace_rho=e_laplace_rhoa, grad_deriv=grad_deriv, & 468 npoints=npoints, epsilon_rho=epsilon_rho, & 469 sx=sx, R=R, gamma=gamma) 470 471 CALL xbecke_roussel_lsd_calc(rho=rhob, norm_drho=norm_drhob, & 472 laplace_rho=laplace_rhob, tau=tau_b, e_0=e_0, e_rho=e_rhob, e_ndrho=e_ndrhob, & 473 e_tau=e_tau_b, e_laplace_rho=e_laplace_rhob, grad_deriv=grad_deriv, & 474 npoints=npoints, epsilon_rho=epsilon_rho, & 475 sx=sx, R=R, gamma=gamma) 476 477!$OMP END PARALLEL 478 479 CALL timestop(handle) 480 END SUBROUTINE xbecke_roussel_lsd_eval 481 482! ************************************************************************************************** 483!> \brief Precalculates which branch of the code has to be taken 484!> There are two main branches, one for a truncated potential and a cutoff 485!> radius, the other for the full coulomb interaction 486!> \param rho grid values 487!> \param norm_drho grid values 488!> \param laplace_rho grid values 489!> \param tau grid values 490!> \param e_0 energies and derivatives 491!> \param e_rho energies and derivatives 492!> \param e_ndrho energies and derivatives 493!> \param e_tau energies and derivatives 494!> \param e_laplace_rho energies and derivatives 495!> \param grad_deriv degree of the derivative that should be evaluated, 496!> if positive all the derivatives up to the given degree are evaluated, 497!> if negative only the given degree is calculated 498!> \param npoints size of the grids 499!> \param epsilon_rho cutoffs 500!> \param sx scales the exchange potential and energies 501!> \param R cutoff Radius for truncated case 502!> \param gamma parameter from original publication, usually set to 1 503!> \par History 504!> 12.2008 created [mguidon] 505!> \author mguidon 506! ************************************************************************************************** 507 SUBROUTINE xbecke_roussel_lsd_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, & 508 e_ndrho, e_tau, e_laplace_rho, grad_deriv, npoints, & 509 epsilon_rho, sx, R, gamma) 510 511 INTEGER, INTENT(in) :: npoints, grad_deriv 512 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_laplace_rho, e_tau, e_ndrho, e_rho, e_0 513 REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: tau, laplace_rho, norm_drho, rho 514 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R, gamma 515 516 CHARACTER(len=*), PARAMETER :: routineN = 'xbecke_roussel_lsd_calc', & 517 routineP = moduleN//':'//routineN 518 519 INTEGER :: ip 520 REAL(dp) :: my_laplace_rho, my_ndrho, my_rho, & 521 my_tau, t1, t15, t16, t2, t3, t4, t5, & 522 t8, t9, yval 523 524! Precalculate y, in order to chose the correct branch afterwards 525 526!$OMP DO 527 528 DO ip = 1, npoints 529 my_rho = MAX(rho(ip), 0.0_dp) 530 IF (my_rho > epsilon_rho) THEN 531 my_ndrho = MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp) 532 my_tau = 1.0_dp*MAX(EPSILON(0.0_dp)*1.e4_dp, tau(ip)) 533 my_laplace_rho = 1.0_dp*laplace_rho(ip) 534 535 t1 = pi**(0.1e1_dp/0.3e1_dp) 536 t2 = t1**2 537 t3 = my_rho**(0.1e1_dp/0.3e1_dp) 538 t4 = t3**2 539 t5 = t4*my_rho 540 t8 = my_ndrho**2 541 t9 = 0.1e1_dp/my_rho 542 ! *** CP2K defines tau in a different way as compared to Becke !!! 543 t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp 544 t16 = 0.1e1_dp/t15 545 yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16 546 IF (R == 0.0_dp) THEN 547 IF (yval <= 0.0_dp) THEN 548 CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), & 549 e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), & 550 sx, gamma, grad_deriv) 551 ELSE 552 CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), & 553 e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), & 554 sx, gamma, grad_deriv) 555 END IF 556 ELSE 557 IF (yval <= 0.0_dp) THEN 558 CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), & 559 e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), & 560 sx, R, gamma, grad_deriv) 561 ELSE 562 CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0(ip), & 563 e_rho(ip), e_ndrho(ip), e_tau(ip), e_laplace_rho(ip), & 564 sx, R, gamma, grad_deriv) 565 END IF 566 END IF 567 END IF 568 END DO 569 570!$OMP END DO 571 572 END SUBROUTINE xbecke_roussel_lsd_calc 573 574! ************************************************************************************************** 575!> \brief full range evaluation for y <= 0 576!> \param rho ... 577!> \param ndrho ... 578!> \param tau ... 579!> \param laplace_rho ... 580!> \param e_0 ... 581!> \param e_rho ... 582!> \param e_ndrho ... 583!> \param e_tau ... 584!> \param e_laplace_rho ... 585!> \param sx ... 586!> \param gamma ... 587!> \param grad_deriv ... 588!> \par History 589!> 12.2008 created [mguidon] 590!> \author mguidon 591! ************************************************************************************************** 592 SUBROUTINE x_br_lsd_y_lte_0(rho, ndrho, tau, laplace_rho, e_0, & 593 e_rho, e_ndrho, e_tau, e_laplace_rho, & 594 sx, gamma, grad_deriv) 595 REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho 596 REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho 597 REAL(dp), INTENT(IN) :: sx, gamma 598 INTEGER, INTENT(IN) :: grad_deriv 599 600 CHARACTER(len=*), PARAMETER :: routineN = 'x_br_lsd_y_lte_0', & 601 routineP = moduleN//':'//routineN 602 603 REAL(KIND=dp) :: t1, t100, t101, t102, t108, t111, t113, t114, t117, t118, t120, t121, t122, & 604 t129, t130, t134, t135, t138, t142, t143, t146, t147, t150, t152, t153, t157, t158, t16, & 605 t161, t164, t165, t166, t169, t17, t170, t172, t173, t19, t199, t2, t20, t202, t207, & 606 t208, t209, t21, t217, t218, t22, t220, t227, t229, t23, t234, t246, t249, t252, t255, & 607 t259, t26, t263, t267, t27, t271, t274, t275, t276, t28, t29, t293, t295, t3, t304, t307, & 608 t308, t32, t320, t33, t34, t340, t341, t342, t344, t346, t349, t35, t350, t353, t354, & 609 t357, t358, t36, t361, t362, t365, t366, t367, t37, t379, t38 610 REAL(KIND=dp) :: t381, t387, t39, t4, t401, t42, t422, t43, t434, t435, t436, t44, t448, & 611 t45, t450, t47, t471, t48, t5, t51, t52, t53, t54, t55, t56, t57, t61, t62, t63, t64, & 612 t66, t67, t70, t71, t72, t75, t78, t81, t84, t87, t88, t89, t9, t91, t92, t93, t94, t95, & 613 t96, t97, yval 614 615 IF (grad_deriv >= 0) THEN 616 t1 = pi**(0.1e1_dp/0.3e1_dp) 617 t2 = t1**2 618 t3 = rho**(0.1e1_dp/0.3e1_dp) 619 t4 = t3**2 620 t5 = t4*rho 621 t9 = ndrho**2 622 t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9/rho/0.4e1_dp)/0.3e1_dp 623 t17 = 0.1e1_dp/t16 624 yval = 0.2e1_dp/0.3e1_dp*t2*t5*t17 625 t19 = t3*rho 626 t20 = 0.3141592654e1_dp**(0.1e1_dp/0.3e1_dp) 627 t21 = t19*t20 628 t22 = br_a1*t2 629 t23 = t5*t17 630 t26 = 0.2e1_dp/0.3e1_dp*t22*t23 + br_a2 631 t27 = ATAN(t26) 632 t28 = -t27 + br_a3 633 t29 = br_c1*t2 634 t32 = t1*pi 635 t33 = br_c2*t32 636 t34 = rho**2 637 t35 = t34*rho 638 t36 = t3*t35 639 t37 = t16**2 640 t38 = 0.1e1_dp/t37 641 t39 = t36*t38 642 t42 = pi**2 643 t43 = br_c3*t42 644 t44 = t34**2 645 t45 = t44*rho 646 t47 = 0.1e1_dp/t37/t16 647 t48 = t45*t47 648 t51 = t2*t42 649 t52 = br_c4*t51 650 t53 = t44*t34 651 t54 = t4*t53 652 t55 = t37**2 653 t56 = 0.1e1_dp/t55 654 t57 = t54*t56 655 t61 = t1*t42*pi 656 t62 = br_c5*t61 657 t63 = t44**2 658 t64 = t3*t63 659 t66 = 0.1e1_dp/t55/t16 660 t67 = t64*t66 661 t70 = br_c0 + 0.2e1_dp/0.3e1_dp*t29*t23 + 0.4e1_dp/0.9e1_dp*t33*t39 & 662 + 0.8e1_dp/0.27e2_dp*t43*t48 + 0.16e2_dp/0.81e2_dp*t52*t57 + 0.32e2_dp & 663 /0.243e3_dp*t62*t67 664 t71 = t28*t70 665 t72 = br_b1*t2 666 t75 = br_b2*t32 667 t78 = br_b3*t42 668 t81 = br_b4*t51 669 t84 = br_b5*t61 670 t87 = br_b0 + 0.2e1_dp/0.3e1_dp*t72*t23 + 0.4e1_dp/0.9e1_dp*t75*t39 & 671 + 0.8e1_dp/0.27e2_dp*t78*t48 + 0.16e2_dp/0.81e2_dp*t81*t57 + 0.32e2_dp & 672 /0.243e3_dp*t84*t67 673 t88 = 0.1e1_dp/t87 674 t89 = t71*t88 675 t91 = EXP(t89/0.3e1_dp) 676 t92 = t21*t91 677 t93 = 0.1e1_dp/t28 678 t94 = 0.1e1_dp/t70 679 t95 = t93*t94 680 t96 = EXP(-t89) 681 t97 = t88*t96 682 t100 = 0.1e1_dp - t96 - t71*t97/0.2e1_dp 683 t101 = t87*t100 684 t102 = t95*t101 685 e_0 = e_0 + (-t92*t102)*sx 686 END IF 687 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 688 t108 = t4*t17 689 t111 = 0.1e1_dp/t3 690 t113 = t38*gamma 691 t114 = t113*t9 692 t117 = 0.10e2_dp/0.9e1_dp*t22*t108 + t22*t111*t114/0.18e2_dp 693 t118 = t26**2 694 t120 = 0.1e1_dp/(0.1e1_dp + t118) 695 t121 = t117*t120 696 t122 = t70*t88 697 t129 = t3*t34 698 t130 = t129*t38 699 t134 = t47*gamma 700 t135 = t134*t9 701 t138 = t44*t47 702 t142 = t56*gamma 703 t143 = t142*t9 704 t146 = t4*t45 705 t147 = t146*t56 706 t150 = t4*t44 707 t152 = t66*gamma 708 t153 = t152*t9 709 t157 = t3*t44*t35 710 t158 = t157*t66 711 t161 = t3*t53 712 t164 = 0.1e1_dp/t55/t37 713 t165 = t164*gamma 714 t166 = t165*t9 715 t169 = 0.10e2_dp/0.9e1_dp*t29*t108 + t29*t111*t114/0.18e2_dp + 0.40e2_dp & 716 /0.27e2_dp*t33*t130 + 0.2e1_dp/0.27e2_dp*t33*t19*t135 + & 717 0.40e2_dp/0.27e2_dp*t43*t138 + 0.2e1_dp/0.27e2_dp*t43*t35*t143 + & 718 0.320e3_dp/0.243e3_dp*t52*t147 + 0.16e2_dp/0.243e3_dp*t52*t150* & 719 t153 + 0.800e3_dp/0.729e3_dp*t62*t158 + 0.40e2_dp/0.729e3_dp*t62*t161 & 720 *t166 721 t170 = t28*t169 722 t172 = t87**2 723 t173 = 0.1e1_dp/t172 724 t199 = 0.10e2_dp/0.9e1_dp*t72*t108 + t72*t111*t114/0.18e2_dp + 0.40e2_dp & 725 /0.27e2_dp*t75*t130 + 0.2e1_dp/0.27e2_dp*t75*t19*t135 + & 726 0.40e2_dp/0.27e2_dp*t78*t138 + 0.2e1_dp/0.27e2_dp*t78*t35*t143 + & 727 0.320e3_dp/0.243e3_dp*t81*t147 + 0.16e2_dp/0.243e3_dp*t81*t150* & 728 t153 + 0.800e3_dp/0.729e3_dp*t84*t158 + 0.40e2_dp/0.729e3_dp*t84*t161 & 729 *t166 730 t202 = -t121*t122 + t170*t88 - t71*t173*t199 731 t207 = t28**2 732 t208 = 0.1e1_dp/t207 733 t209 = t91*t208 734 t217 = t21*t91*t93 735 t218 = t70**2 736 t220 = 0.1e1_dp/t218*t87 737 t227 = -t202 738 t229 = t122*t96 739 t234 = t173*t96 740 e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t3*t20*t91*t102 - t21*t202*t91* & 741 t102/0.3e1_dp - t21*t209*t94*t87*t100*t117*t120 + t217 & 742 *t220*t100*t169 - t92*t95*t199*t100 - t92*t95*t87* & 743 (-t227*t96 + t121*t229/0.2e1_dp - t170*t97/0.2e1_dp + t71*t234 & 744 *t199/0.2e1_dp - t71*t88*t227*t96/0.2e1_dp))*sx 745 t246 = t4*t38 746 t249 = t120*t70 747 t252 = t22*t246*gamma*ndrho*t249*t88 748 t255 = t113*ndrho 749 t259 = t134*ndrho 750 t263 = t142*ndrho 751 t267 = t152*ndrho 752 t271 = t165*ndrho 753 t274 = -t29*t4*t255/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t33*t129*t259 & 754 - 0.4e1_dp/0.27e2_dp*t43*t44*t263 - 0.32e2_dp/0.243e3_dp*t52*t146 & 755 *t267 - 0.80e2_dp/0.729e3_dp*t62*t157*t271 756 t275 = t28*t274 757 t276 = t275*t88 758 t293 = -t72*t4*t255/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t75*t129*t259 & 759 - 0.4e1_dp/0.27e2_dp*t78*t44*t263 - 0.32e2_dp/0.243e3_dp*t81*t146 & 760 *t267 - 0.80e2_dp/0.729e3_dp*t84*t157*t271 761 t295 = t71*t173*t293 762 t304 = t208*t94*t87 763 t307 = t100*br_a1*t2 764 t308 = ndrho*t120 765 t320 = -t252/0.9e1_dp - t276 + t295 766 e_ndrho = e_ndrho + (-t21*(t252/0.27e2_dp + t276/0.3e1_dp - t295/0.3e1_dp)*t91 & 767 *t102 + t34*t20*t91*t304*t307*t113*t308/0.9e1_dp + t217 & 768 *t220*t100*t274 - t92*t95*t293*t100 - t92*t95*t87 & 769 *(-t320*t96 - t22*t246*gamma*t308*t229/0.18e2_dp - t275 & 770 *t97/0.2e1_dp + t71*t234*t293/0.2e1_dp - t71*t88*t320*t96 & 771 /0.2e1_dp))*sx 772 t340 = t5*t38 773 t341 = t22*t340 774 t342 = gamma*t120 775 t344 = t341*t342*t122 776 t346 = t340*gamma 777 t349 = t36*t47 778 t350 = t349*gamma 779 t353 = t45*t56 780 t354 = t353*gamma 781 t357 = t54*t66 782 t358 = t357*gamma 783 t361 = t64*t164 784 t362 = t361*gamma 785 t365 = 0.4e1_dp/0.9e1_dp*t29*t346 + 0.16e2_dp/0.27e2_dp*t33*t350 + & 786 0.16e2_dp/0.27e2_dp*t43*t354 + 0.128e3_dp/0.243e3_dp*t52*t358 + 0.320e3_dp & 787 /0.729e3_dp*t62*t362 788 t366 = t28*t365 789 t367 = t366*t88 790 t379 = 0.4e1_dp/0.9e1_dp*t72*t346 + 0.16e2_dp/0.27e2_dp*t75*t350 + & 791 0.16e2_dp/0.27e2_dp*t78*t354 + 0.128e3_dp/0.243e3_dp*t81*t358 + 0.320e3_dp & 792 /0.729e3_dp*t84*t362 793 t381 = t71*t173*t379 794 t387 = t35*t20 795 t401 = 0.4e1_dp/0.9e1_dp*t344 - t367 + t381 796 e_tau = e_tau + (-t21*(-0.4e1_dp/0.27e2_dp*t344 + t367/0.3e1_dp - t381/0.3e1_dp) & 797 *t91*t102 - 0.4e1_dp/0.9e1_dp*t387*t91*t304*t307*t113* & 798 t120 + t217*t220*t100*t365 - t92*t95*t379*t100 - t92 & 799 *t95*t87*(-t401*t96 + 0.2e1_dp/0.9e1_dp*t341*t342*t229 - & 800 t366*t97/0.2e1_dp + t71*t234*t379/0.2e1_dp - t71*t88*t401 & 801 *t96/0.2e1_dp))*sx 802 t422 = t22*t5*t38*t120*t122 803 t434 = -t29*t340/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t33*t349 - 0.4e1_dp/ & 804 0.27e2_dp*t43*t353 - 0.32e2_dp/0.243e3_dp*t52*t357 - 0.80e2_dp/0.729e3_dp & 805 *t62*t361 806 t435 = t28*t434 807 t436 = t435*t88 808 t448 = -t72*t340/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t75*t349 - 0.4e1_dp/ & 809 0.27e2_dp*t78*t353 - 0.32e2_dp/0.243e3_dp*t81*t357 - 0.80e2_dp/0.729e3_dp & 810 *t84*t361 811 t450 = t71*t173*t448 812 t471 = -t422/0.9e1_dp - t436 + t450 813 e_laplace_rho = e_laplace_rho + (-t21*(t422/0.27e2_dp + t436/0.3e1_dp - t450/0.3e1_dp)*t91* & 814 t102 + t387*t209*t94*t101*br_a1*t2*t38*t120/0.9e1_dp & 815 + t217*t220*t100*t434 - t92*t95*t448*t100 - t92*t95 & 816 *t87*(-t471*t96 - t341*t249*t97/0.18e2_dp - t435*t97/0.2e1_dp & 817 + t71*t234*t448/0.2e1_dp - t71*t88*t471*t96/0.2e1_dp))*sx 818 END IF 819 END SUBROUTINE x_br_lsd_y_lte_0 820 821! ************************************************************************************************** 822!> \brief Full range evaluation for y > 0 823!> \param rho ... 824!> \param ndrho ... 825!> \param tau ... 826!> \param laplace_rho ... 827!> \param e_0 ... 828!> \param e_rho ... 829!> \param e_ndrho ... 830!> \param e_tau ... 831!> \param e_laplace_rho ... 832!> \param sx ... 833!> \param gamma ... 834!> \param grad_deriv ... 835!> \par History 836!> 12.2008 created [mguidon] 837!> \author mguidon 838! ************************************************************************************************** 839 SUBROUTINE x_br_lsd_y_gt_0(rho, ndrho, tau, laplace_rho, e_0, & 840 e_rho, e_ndrho, e_tau, e_laplace_rho, & 841 sx, gamma, grad_deriv) 842 REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho 843 REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho 844 REAL(dp), INTENT(IN) :: sx, gamma 845 INTEGER, INTENT(IN) :: grad_deriv 846 847 CHARACTER(len=*), PARAMETER :: routineN = 'x_br_lsd_y_gt_0', & 848 routineP = moduleN//':'//routineN 849 850 REAL(KIND=dp) :: t1, t102, t103, t104, t106, t107, t108, t109, t110, t111, t112, t115, t117, & 851 t124, t131, t134, t137, t138, t142, t143, t154, t157, t158, t159, t16, t160, t162, t165, & 852 t167, t168, t17, t176, t180, t181, t184, t185, t188, t19, t190, t191, t195, t196, t199, & 853 t2, t20, t202, t203, t204, t207, t208, t21, t210, t211, t22, t23, t237, t24, t240, t245, & 854 t248, t249, t25, t255, t256, t258, t26, t265, t267, t272, t285, t288, t289, t29, t297, & 855 t298, t3, t30, t301, t305, t309, t31, t313, t317, t32, t320, t321, t33, t338, t34, t341, & 856 t35, t356, t36, t37, t376, t377, t382, t383, t387, t388, t391 857 REAL(KIND=dp) :: t392, t395, t396, t399, t4, t400, t403, t404, t41, t416, t419, t42, t43, & 858 t434, t458, t459, t47, t471, t472, t48, t484, t487, t49, t5, t50, t502, t51, t54, t57, & 859 t58, t59, t6, t60, t62, t63, t66, t67, t68, t69, t70, t71, t72, t76, t77, t78, t79, t81, & 860 t82, t85, t86, t87, t9, t90, t93, t96, t99, yval 861 862 IF (grad_deriv >= 0) THEN 863 t1 = pi**(0.1e1_dp/0.3e1_dp) 864 t2 = t1**2 865 t3 = rho**(0.1e1_dp/0.3e1_dp) 866 t4 = t3**2 867 t5 = t4*rho 868 t6 = t2*t5 869 t9 = ndrho**2 870 t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9/rho/0.4e1_dp)/0.3e1_dp 871 t17 = 0.1e1_dp/t16 872 yval = 0.2e1_dp/0.3e1_dp*t6*t17 873 t19 = t3*rho 874 t20 = 0.3141592654e1_dp**(0.1e1_dp/0.3e1_dp) 875 t21 = t19*t20 876 t22 = 0.1e1_dp/br_BB 877 t23 = 0.1e1_dp/t2 878 t24 = t22*t23 879 t25 = 0.1e1_dp/t5 880 t26 = t25*t16 881 t29 = br_BB**2 882 t30 = t1*pi 883 t31 = t29*t30 884 t32 = rho**2 885 t33 = t32*rho 886 t34 = t3*t33 887 t35 = t16**2 888 t36 = 0.1e1_dp/t35 889 t37 = t34*t36 890 t41 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t31*t37) 891 t42 = t41*t22 892 t43 = t23*t25 893 t47 = 0.1500000000000000e1_dp*t24*t26 + 0.3e1_dp/0.2e1_dp*t42*t43 & 894 *t16 895 t48 = LOG(t47) 896 t49 = t48 + 0.2e1_dp 897 t50 = br_d1*t2 898 t51 = t5*t17 899 t54 = br_d2*t30 900 t57 = pi**2 901 t58 = br_d3*t57 902 t59 = t32**2 903 t60 = t59*rho 904 t62 = 0.1e1_dp/t35/t16 905 t63 = t60*t62 906 t66 = t2*t57 907 t67 = br_d4*t66 908 t68 = t59*t32 909 t69 = t4*t68 910 t70 = t35**2 911 t71 = 0.1e1_dp/t70 912 t72 = t69*t71 913 t76 = t1*t57*pi 914 t77 = br_d5*t76 915 t78 = t59**2 916 t79 = t3*t78 917 t81 = 0.1e1_dp/t70/t16 918 t82 = t79*t81 919 t85 = br_d0 + 0.2e1_dp/0.3e1_dp*t50*t51 + 0.4e1_dp/0.9e1_dp*t54*t37 & 920 + 0.8e1_dp/0.27e2_dp*t58*t63 + 0.16e2_dp/0.81e2_dp*t67*t72 + 0.32e2_dp & 921 /0.243e3_dp*t77*t82 922 t86 = t49*t85 923 t87 = br_e1*t2 924 t90 = br_e2*t30 925 t93 = br_e3*t57 926 t96 = br_e4*t66 927 t99 = br_e5*t76 928 t102 = br_e0 + 0.2e1_dp/0.3e1_dp*t87*t51 + 0.4e1_dp/0.9e1_dp*t90*t37 & 929 + 0.8e1_dp/0.27e2_dp*t93*t63 + 0.16e2_dp/0.81e2_dp*t96*t72 + 0.32e2_dp & 930 /0.243e3_dp*t99*t82 931 t103 = 0.1e1_dp/t102 932 t104 = t86*t103 933 t106 = EXP(t104/0.3e1_dp) 934 t107 = t21*t106 935 t108 = 0.1e1_dp/t49 936 t109 = 0.1e1_dp/t85 937 t110 = t108*t109 938 t111 = EXP(-t104) 939 t112 = t103*t111 940 t115 = 0.1e1_dp - t111 - t86*t112/0.2e1_dp 941 t117 = t110*t102*t115 942 e_0 = e_0 + (-t107*t117)*sx 943 END IF 944 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 945 t124 = 0.1e1_dp/t4/t32 946 t131 = 0.1e1_dp/t4/t33*gamma*t9 947 t134 = 0.1e1_dp/t41 948 t137 = t3*t32 949 t138 = t137*t36 950 t142 = t62*gamma 951 t143 = t142*t9 952 t154 = t42*t23 953 t157 = -0.2500000000e1_dp*t24*t124*t16 - 0.1250000000e0_dp*t24* & 954 t131 + 0.3e1_dp/0.4e1_dp*t134*t22*t23*t26*(0.40e2_dp/0.27e2_dp* & 955 t31*t138 + 0.2e1_dp/0.27e2_dp*t31*t19*t143) - 0.5e1_dp/0.2e1_dp* & 956 t42*t23*t124*t16 - t154*t131/0.8e1_dp 957 t158 = 0.1e1_dp/t47 958 t159 = t157*t158 959 t160 = t85*t103 960 t162 = t4*t17 961 t165 = 0.1e1_dp/t3 962 t167 = t36*gamma 963 t168 = t167*t9 964 t176 = t59*t62 965 t180 = t71*gamma 966 t181 = t180*t9 967 t184 = t4*t60 968 t185 = t184*t71 969 t188 = t4*t59 970 t190 = t81*gamma 971 t191 = t190*t9 972 t195 = t3*t59*t33 973 t196 = t195*t81 974 t199 = t3*t68 975 t202 = 0.1e1_dp/t70/t35 976 t203 = t202*gamma 977 t204 = t203*t9 978 t207 = 0.10e2_dp/0.9e1_dp*t50*t162 + t50*t165*t168/0.18e2_dp + 0.40e2_dp & 979 /0.27e2_dp*t54*t138 + 0.2e1_dp/0.27e2_dp*t54*t19*t143 + & 980 0.40e2_dp/0.27e2_dp*t58*t176 + 0.2e1_dp/0.27e2_dp*t58*t33*t181 + & 981 0.320e3_dp/0.243e3_dp*t67*t185 + 0.16e2_dp/0.243e3_dp*t67*t188* & 982 t191 + 0.800e3_dp/0.729e3_dp*t77*t196 + 0.40e2_dp/0.729e3_dp*t77*t199 & 983 *t204 984 t208 = t49*t207 985 t210 = t102**2 986 t211 = 0.1e1_dp/t210 987 t237 = 0.10e2_dp/0.9e1_dp*t87*t162 + t87*t165*t168/0.18e2_dp + 0.40e2_dp & 988 /0.27e2_dp*t90*t138 + 0.2e1_dp/0.27e2_dp*t90*t19*t143 + & 989 0.40e2_dp/0.27e2_dp*t93*t176 + 0.2e1_dp/0.27e2_dp*t93*t33*t181 + & 990 0.320e3_dp/0.243e3_dp*t96*t185 + 0.16e2_dp/0.243e3_dp*t96*t188* & 991 t191 + 0.800e3_dp/0.729e3_dp*t99*t196 + 0.40e2_dp/0.729e3_dp*t99*t199 & 992 *t204 993 t240 = t159*t160 + t208*t103 - t86*t211*t237 994 t245 = t49**2 995 t248 = t21*t106/t245 996 t249 = t109*t102 997 t255 = t21*t106*t108 998 t256 = t85**2 999 t258 = 0.1e1_dp/t256*t102 1000 t265 = -t240 1001 t267 = t160*t111 1002 t272 = t211*t111 1003 e_rho = e_rho + (-0.4e1_dp/0.3e1_dp*t3*t20*t106*t117 - t21*t240*t106 & 1004 *t117/0.3e1_dp + t248*t249*t115*t157*t158 + t255*t258* & 1005 t115*t207 - t107*t110*t237*t115 - t107*t110*t102*(-t265 & 1006 *t111 - t159*t267/0.2e1_dp - t208*t112/0.2e1_dp + t86*t272 & 1007 *t237/0.2e1_dp - t86*t103*t265*t111/0.2e1_dp))*sx 1008 t285 = t124*gamma*ndrho 1009 t288 = t134*br_BB 1010 t289 = t288*t2 1011 t297 = 0.2500000000000000e0_dp*t24*t285 - t289*t4*t36*gamma* & 1012 ndrho/0.9e1_dp + t154*t285/0.4e1_dp 1013 t298 = t297*t158 1014 t301 = t167*ndrho 1015 t305 = t142*ndrho 1016 t309 = t180*ndrho 1017 t313 = t190*ndrho 1018 t317 = t203*ndrho 1019 t320 = -t50*t4*t301/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t54*t137*t305 & 1020 - 0.4e1_dp/0.27e2_dp*t58*t59*t309 - 0.32e2_dp/0.243e3_dp*t67*t184 & 1021 *t313 - 0.80e2_dp/0.729e3_dp*t77*t195*t317 1022 t321 = t49*t320 1023 t338 = -t87*t4*t301/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t90*t137*t305 & 1024 - 0.4e1_dp/0.27e2_dp*t93*t59*t309 - 0.32e2_dp/0.243e3_dp*t96*t184 & 1025 *t313 - 0.80e2_dp/0.729e3_dp*t99*t195*t317 1026 t341 = t298*t160 + t321*t103 - t86*t211*t338 1027 t356 = -t341 1028 e_ndrho = e_ndrho + (-t21*t341*t106*t117/0.3e1_dp + t248*t249*t115*t297 & 1029 *t158 + t255*t258*t115*t320 - t107*t110*t338*t115 & 1030 - t107*t110*t102*(-t356*t111 - t298*t267/0.2e1_dp - t321 & 1031 *t112/0.2e1_dp + t86*t272*t338/0.2e1_dp - t86*t103*t356* & 1032 t111/0.2e1_dp))*sx 1033 t376 = t5*t36 1034 t377 = t376*gamma 1035 t382 = -0.1000000000e1_dp*t24*t25*gamma + 0.4e1_dp/0.9e1_dp*t289*t377 & 1036 - t42*t43*gamma 1037 t383 = t382*t158 1038 t387 = t34*t62 1039 t388 = t387*gamma 1040 t391 = t60*t71 1041 t392 = t391*gamma 1042 t395 = t69*t81 1043 t396 = t395*gamma 1044 t399 = t79*t202 1045 t400 = t399*gamma 1046 t403 = 0.4e1_dp/0.9e1_dp*t50*t377 + 0.16e2_dp/0.27e2_dp*t54*t388 + & 1047 0.16e2_dp/0.27e2_dp*t58*t392 + 0.128e3_dp/0.243e3_dp*t67*t396 + 0.320e3_dp & 1048 /0.729e3_dp*t77*t400 1049 t404 = t49*t403 1050 t416 = 0.4e1_dp/0.9e1_dp*t87*t377 + 0.16e2_dp/0.27e2_dp*t90*t388 + & 1051 0.16e2_dp/0.27e2_dp*t93*t392 + 0.128e3_dp/0.243e3_dp*t96*t396 + 0.320e3_dp & 1052 /0.729e3_dp*t99*t400 1053 t419 = t383*t160 + t404*t103 - t86*t211*t416 1054 t434 = -t419 1055 e_tau = e_tau + (-t21*t419*t106*t117/0.3e1_dp + t248*t249*t115*t382 & 1056 *t158 + t255*t258*t115*t403 - t107*t110*t416*t115 - & 1057 t107*t110*t102*(-t434*t111 - t383*t267/0.2e1_dp - t404* & 1058 t112/0.2e1_dp + t86*t272*t416/0.2e1_dp - t86*t103*t434*t111 & 1059 /0.2e1_dp))*sx 1060 t458 = 0.2500000000000000e0_dp*t24*t25 - t288*t6*t36/0.9e1_dp + & 1061 t42*t43/0.4e1_dp 1062 t459 = t458*t158 1063 t471 = -t50*t376/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t54*t387 - 0.4e1_dp/ & 1064 0.27e2_dp*t58*t391 - 0.32e2_dp/0.243e3_dp*t67*t395 - 0.80e2_dp/0.729e3_dp & 1065 *t77*t399 1066 t472 = t49*t471 1067 t484 = -t87*t376/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t90*t387 - 0.4e1_dp/ & 1068 0.27e2_dp*t93*t391 - 0.32e2_dp/0.243e3_dp*t96*t395 - 0.80e2_dp/0.729e3_dp & 1069 *t99*t399 1070 t487 = t459*t160 + t472*t103 - t86*t211*t484 1071 t502 = -t487 1072 e_laplace_rho = e_laplace_rho + (-t21*t487*t106*t117/0.3e1_dp + t248*t249*t115*t458 & 1073 *t158 + t255*t258*t115*t471 - t107*t110*t484*t115 - & 1074 t107*t110*t102*(-t502*t111 - t459*t267/0.2e1_dp - t472* & 1075 t112/0.2e1_dp + t86*t272*t484/0.2e1_dp - t86*t103*t502*t111 & 1076 /0.2e1_dp))*sx 1077 END IF 1078 1079 END SUBROUTINE x_br_lsd_y_gt_0 1080 1081! ************************************************************************************************** 1082!> \brief Truncated long range part for y <= 0 1083!> \param rho ... 1084!> \param ndrho ... 1085!> \param tau ... 1086!> \param laplace_rho ... 1087!> \param e_0 ... 1088!> \param e_rho ... 1089!> \param e_ndrho ... 1090!> \param e_tau ... 1091!> \param e_laplace_rho ... 1092!> \param sx ... 1093!> \param R ... 1094!> \param gamma ... 1095!> \param grad_deriv ... 1096!> \par History 1097!> 12.2008 created [mguidon] 1098!> \author mguidon 1099! ************************************************************************************************** 1100 SUBROUTINE x_br_lsd_y_lte_0_cutoff(rho, ndrho, tau, laplace_rho, e_0, & 1101 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1102 sx, R, gamma, grad_deriv) 1103 REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho 1104 REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho 1105 REAL(dp), INTENT(IN) :: sx, R, gamma 1106 INTEGER, INTENT(IN) :: grad_deriv 1107 1108 CHARACTER(len=*), PARAMETER :: routineN = 'x_br_lsd_y_lte_0_cutoff', & 1109 routineP = moduleN//':'//routineN 1110 1111 REAL(KIND=dp) :: brval, t1, t101, t11, t12, t18, t2, t20, & 1112 t24, t25, t26, t3, t31, t33, t36, t38, & 1113 t4, t41, t43, t47, t50, t54, t56, t6, & 1114 t60, t62, t66, t69, t7, t70, t88, t89, & 1115 t96 1116 1117 t1 = 8**(0.1e1_dp/0.3e1_dp) 1118 t2 = t1**2 1119 t3 = pi**(0.1e1_dp/0.3e1_dp) 1120 t4 = t3**2 1121 t6 = rho**(0.1e1_dp/0.3e1_dp) 1122 t7 = t6**2 1123 t11 = ndrho**2 1124 t12 = 0.1e1_dp/rho 1125 t18 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t11*t12/0.4e1_dp)/0.3e1_dp 1126 t20 = t7*rho/t18 1127 t24 = ATAN(0.2e1_dp/0.3e1_dp*br_a1*t4*t20 + br_a2) 1128 t25 = -t24 + br_a3 1129 t26 = t25**2 1130 t31 = t3*pi 1131 t33 = rho**2 1132 t36 = t18**2 1133 t38 = t6*t33*rho/t36 1134 t41 = pi**2 1135 t43 = t33**2 1136 t47 = t43*rho/t36/t18 1137 t50 = t4*t41 1138 t54 = t36**2 1139 t56 = t7*t43*t33/t54 1140 t60 = t3*t41*pi 1141 t62 = t43**2 1142 t66 = t6*t62/t54/t18 1143 t69 = br_c0 + 0.2e1_dp/0.3e1_dp*br_c1*t4*t20 + 0.4e1_dp/0.9e1_dp*br_c2 & 1144 *t31*t38 + 0.8e1_dp/0.27e2_dp*br_c3*t41*t47 + 0.16e2_dp/0.81e2_dp & 1145 *br_c4*t50*t56 + 0.32e2_dp/0.243e3_dp*br_c5*t60*t66 1146 t70 = t69**2 1147 t88 = br_b0 + 0.2e1_dp/0.3e1_dp*br_b1*t4*t20 + 0.4e1_dp/0.9e1_dp*br_b2 & 1148 *t31*t38 + 0.8e1_dp/0.27e2_dp*br_b3*t41*t47 + 0.16e2_dp/0.81e2_dp & 1149 *br_b4*t50*t56 + 0.32e2_dp/0.243e3_dp*br_b5*t60*t66 1150 t89 = t88**2 1151 t96 = EXP(-t25*t69/t88) 1152 t101 = (t26*t25*t70*t69/t89/t88*t96/0.3141592654e1_dp* & 1153 t12)**(0.1e1_dp/0.3e1_dp) 1154 brval = REAL(t2, KIND=dp)*t101/0.8e1_dp 1155 1156 IF (R > brval) THEN 1157 CALL x_br_lsd_y_lte_0_cutoff_R_gt_b(rho, ndrho, tau, laplace_rho, e_0, & 1158 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1159 sx, R, gamma, grad_deriv) 1160 ELSE 1161 CALL x_br_lsd_y_lte_0_cutoff_R_lte_b(rho, ndrho, tau, laplace_rho, e_0, & 1162 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1163 sx, R, gamma, grad_deriv) 1164 END IF 1165 1166 END SUBROUTINE x_br_lsd_y_lte_0_cutoff 1167 1168! ************************************************************************************************** 1169!> \brief Truncated long range part for y <= 0 and R > b 1170!> \param rho ... 1171!> \param ndrho ... 1172!> \param tau ... 1173!> \param laplace_rho ... 1174!> \param e_0 ... 1175!> \param e_rho ... 1176!> \param e_ndrho ... 1177!> \param e_tau ... 1178!> \param e_laplace_rho ... 1179!> \param sx ... 1180!> \param R ... 1181!> \param gamma ... 1182!> \param grad_deriv ... 1183!> \par History 1184!> 12.2008 created [mguidon] 1185!> \author mguidon 1186! ************************************************************************************************** 1187 SUBROUTINE x_br_lsd_y_lte_0_cutoff_R_gt_b(rho, ndrho, tau, laplace_rho, e_0, & 1188 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1189 sx, R, gamma, grad_deriv) 1190 REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho 1191 REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho 1192 REAL(dp), INTENT(IN) :: sx, R, gamma 1193 INTEGER, INTENT(IN) :: grad_deriv 1194 1195 CHARACTER(len=*), PARAMETER :: routineN = 'x_br_lsd_y_lte_0_cutoff_R_gt_b', & 1196 routineP = moduleN//':'//routineN 1197 1198 REAL(KIND=dp) :: t1, t10, t100, t101, t102, t103, t104, t106, t108, t109, t110, t114, t117, & 1199 t119, t120, t123, t124, t129, t132, t134, t135, t138, t139, t141, t142, t143, t149, t150, & 1200 t153, t155, t156, t159, t16, t163, t164, t167, t168, t17, t171, t173, t174, t178, t179, & 1201 t18, t182, t185, t186, t187, t190, t192, t193, t2, t21, t219, t22, t221, t223, t224, & 1202 t225, t227, t228, t23, t231, t232, t233, t235, t24, t240, t245, t247, t259, t261, t263, & 1203 t265, t268, t269, t27, t270, t272, t274, t275, t28, t283, t288, t29, t291, t3, t30, t301, & 1204 t305, t31, t312, t314, t315, t317, t319, t32, t323, t327, t33 1205 REAL(KIND=dp) :: t331, t335, t338, t34, t340, t356, t358, t359, t361, t363, t364, t366, & 1206 t367, t37, t38, t388, t39, t390, t392, t394, t397, t399, t4, t40, t403, t405, t409, t410, & 1207 t414, t42, t423, t426, t43, t434, t443, t450, t455, t456, t459, t46, t460, t463, t464, & 1208 t467, t468, t47, t471, t472, t475, t477, t48, t488, t49, t490, t493, t494, t496, t497, & 1209 t499, t5, t50, t51, t516, t518, t52, t520, t522, t525, t526, t527, t530, t532, t545, & 1210 t548, t56, t563, t57, t572, t574, t58, t585, t587, t59, t598, t6, t600, t605, t606, t608, & 1211 t609, t61, t62, t628, t630, t632, t634, t639, t641, t645, t65, t655 1212 REAL(KIND=dp) :: t66, t67, t672, t70, t73, t76, t79, t82, t83, t84, t85, t86, t87, t88, t89, & 1213 t9, t90, t91, t93, t94, t95, t96, t97, t99 1214 1215 IF (grad_deriv >= 0) THEN 1216 t1 = pi**(0.1e1_dp/0.3e1_dp) 1217 t2 = t1**2 1218 t3 = br_a1*t2 1219 t4 = rho**(0.1e1_dp/0.3e1_dp) 1220 t5 = t4**2 1221 t6 = t5*rho 1222 t9 = ndrho**2 1223 t10 = 0.1e1_dp/rho 1224 t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9*t10/0.4e1_dp)/0.3e1_dp 1225 t17 = 0.1e1_dp/t16 1226 t18 = t6*t17 1227 t21 = 0.2e1_dp/0.3e1_dp*t3*t18 + br_a2 1228 t22 = ATAN(t21) 1229 t23 = -t22 + br_a3 1230 t24 = br_c1*t2 1231 t27 = t1*pi 1232 t28 = br_c2*t27 1233 t29 = rho**2 1234 t30 = t29*rho 1235 t31 = t4*t30 1236 t32 = t16**2 1237 t33 = 0.1e1_dp/t32 1238 t34 = t31*t33 1239 t37 = pi**2 1240 t38 = br_c3*t37 1241 t39 = t29**2 1242 t40 = t39*rho 1243 t42 = 0.1e1_dp/t32/t16 1244 t43 = t40*t42 1245 t46 = t2*t37 1246 t47 = br_c4*t46 1247 t48 = t39*t29 1248 t49 = t5*t48 1249 t50 = t32**2 1250 t51 = 0.1e1_dp/t50 1251 t52 = t49*t51 1252 t56 = t1*t37*pi 1253 t57 = br_c5*t56 1254 t58 = t39**2 1255 t59 = t4*t58 1256 t61 = 0.1e1_dp/t50/t16 1257 t62 = t59*t61 1258 t65 = br_c0 + 0.2e1_dp/0.3e1_dp*t24*t18 + 0.4e1_dp/0.9e1_dp*t28*t34 & 1259 + 0.8e1_dp/0.27e2_dp*t38*t43 + 0.16e2_dp/0.81e2_dp*t47*t52 + 0.32e2_dp & 1260 /0.243e3_dp*t57*t62 1261 t66 = t23*t65 1262 t67 = br_b1*t2 1263 t70 = br_b2*t27 1264 t73 = br_b3*t37 1265 t76 = br_b4*t46 1266 t79 = br_b5*t56 1267 t82 = br_b0 + 0.2e1_dp/0.3e1_dp*t67*t18 + 0.4e1_dp/0.9e1_dp*t70*t34 & 1268 + 0.8e1_dp/0.27e2_dp*t73*t43 + 0.16e2_dp/0.81e2_dp*t76*t52 + 0.32e2_dp & 1269 /0.243e3_dp*t79*t62 1270 t83 = 0.1e1_dp/t82 1271 t84 = t66*t83 1272 t85 = 8**(0.1e1_dp/0.3e1_dp) 1273 t86 = t23**2 1274 t87 = t86*t23 1275 t88 = t65**2 1276 t89 = t88*t65 1277 t90 = t87*t89 1278 t91 = t82**2 1279 t93 = 0.1e1_dp/t91/t82 1280 t94 = t90*t93 1281 t95 = EXP(-t84) 1282 t96 = 0.1e1_dp/0.3141592654e1_dp 1283 t97 = t95*t96 1284 t99 = t94*t97*t10 1285 t100 = t99**(0.1e1_dp/0.3e1_dp) 1286 t101 = 0.1e1_dp/t100 1287 t102 = REAL(t85, KIND=dp)*t101 1288 t103 = t102*R 1289 t104 = t84*t103 1290 t106 = EXP(t84 - t104) 1291 t108 = t106*t23 1292 t109 = t65*t83 1293 t110 = t108*t109 1294 t114 = t83*REAL(t85, KIND=dp)*t101*R 1295 t117 = EXP(-t84 - t104) 1296 t119 = t117*t23 1297 t120 = t119*t109 1298 t123 = -0.2e1_dp*t106 + t110 - t108*t65*t114 + 0.2e1_dp*t117 + t120 & 1299 + t119*t65*t114 1300 t124 = rho*t123 1301 e_0 = e_0 + (t124*t102/0.8e1_dp)*sx 1302 END IF 1303 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 1304 t129 = t5*t17 1305 t132 = 0.1e1_dp/t4 1306 t134 = t33*gamma 1307 t135 = t134*t9 1308 t138 = 0.10e2_dp/0.9e1_dp*t3*t129 + t3*t132*t135/0.18e2_dp 1309 t139 = t21**2 1310 t141 = 0.1e1_dp/(0.1e1_dp + t139) 1311 t142 = t138*t141 1312 t143 = t142*t109 1313 t149 = t4*t29 1314 t150 = t149*t33 1315 t153 = t4*rho 1316 t155 = t42*gamma 1317 t156 = t155*t9 1318 t159 = t39*t42 1319 t163 = t51*gamma 1320 t164 = t163*t9 1321 t167 = t5*t40 1322 t168 = t167*t51 1323 t171 = t5*t39 1324 t173 = t61*gamma 1325 t174 = t173*t9 1326 t178 = t4*t39*t30 1327 t179 = t178*t61 1328 t182 = t4*t48 1329 t185 = 0.1e1_dp/t50/t32 1330 t186 = t185*gamma 1331 t187 = t186*t9 1332 t190 = 0.10e2_dp/0.9e1_dp*t24*t129 + t24*t132*t135/0.18e2_dp + 0.40e2_dp & 1333 /0.27e2_dp*t28*t150 + 0.2e1_dp/0.27e2_dp*t28*t153*t156 + & 1334 0.40e2_dp/0.27e2_dp*t38*t159 + 0.2e1_dp/0.27e2_dp*t38*t30*t164 & 1335 + 0.320e3_dp/0.243e3_dp*t47*t168 + 0.16e2_dp/0.243e3_dp*t47*t171* & 1336 t174 + 0.800e3_dp/0.729e3_dp*t57*t179 + 0.40e2_dp/0.729e3_dp*t57* & 1337 t182*t187 1338 t192 = t23*t190*t83 1339 t193 = 0.1e1_dp/t91 1340 t219 = 0.10e2_dp/0.9e1_dp*t67*t129 + t67*t132*t135/0.18e2_dp + 0.40e2_dp & 1341 /0.27e2_dp*t70*t150 + 0.2e1_dp/0.27e2_dp*t70*t153*t156 + & 1342 0.40e2_dp/0.27e2_dp*t73*t159 + 0.2e1_dp/0.27e2_dp*t73*t30*t164 & 1343 + 0.320e3_dp/0.243e3_dp*t76*t168 + 0.16e2_dp/0.243e3_dp*t76*t171* & 1344 t174 + 0.800e3_dp/0.729e3_dp*t79*t179 + 0.40e2_dp/0.729e3_dp*t79* & 1345 t182*t187 1346 t221 = t66*t193*t219 1347 t223 = t142*t65*t114 1348 t224 = t192*t103 1349 t225 = t66*t193 1350 t227 = t102*R*t219 1351 t228 = t225*t227 1352 t231 = REAL(t85, KIND=dp)/t100/t99 1353 t232 = t86*t89 1354 t233 = t93*t95 1355 t235 = t96*t10 1356 t240 = t87*t88*t93 1357 t245 = t91**2 1358 t247 = t90/t245 1359 t259 = -0.3e1_dp*t232*t233*t235*t142 + 0.3e1_dp*t240*t97*t10 & 1360 *t190 - 0.3e1_dp*t247*t97*t10*t219 + t94*(t143 - t192 + & 1361 t221)*t95*t235 - t94*t97/t29 1362 t261 = t231*R*t259 1363 t263 = t84*t261/0.3e1_dp 1364 t265 = (-t143 + t192 - t221 + t223 - t224 + t228 + t263)*t106 1365 t268 = t106*t138 1366 t269 = t141*t65 1367 t270 = t269*t83 1368 t272 = t190*t83 1369 t274 = t65*t193 1370 t275 = t274*t219 1371 t283 = t108*t274 1372 t288 = (t143 - t192 + t221 + t223 - t224 + t228 + t263)*t117 1373 t291 = t117*t138 1374 t301 = t119*t274 1375 t305 = -0.2e1_dp*t265 + t265*t84 - t268*t270 + t108*t272 - t108 & 1376 *t275 - t265*t66*t114 + t268*t269*t114 - t108*t190* & 1377 t114 + t283*t227 + t110*t261/0.3e1_dp + 0.2e1_dp*t288 + t288*t84 & 1378 - t291*t270 + t119*t272 - t119*t275 + t288*t66*t114 - & 1379 t291*t269*t114 + t119*t190*t114 - t301*t227 - t120*t261 & 1380 /0.3e1_dp 1381 e_rho = e_rho + (t123*REAL(t85, KIND=dp)*t101/0.8e1_dp + rho*t305*t102/0.8e1_dp & 1382 - t124*t231*t259/0.24e2_dp)*sx 1383 t312 = t5*t33 1384 t314 = gamma*ndrho 1385 t315 = t314*t270 1386 t317 = t3*t312*t315/0.9e1_dp 1387 t319 = t134*ndrho 1388 t323 = t155*ndrho 1389 t327 = t163*ndrho 1390 t331 = t173*ndrho 1391 t335 = t186*ndrho 1392 t338 = -t24*t5*t319/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t28*t149*t323 & 1393 - 0.4e1_dp/0.27e2_dp*t38*t39*t327 - 0.32e2_dp/0.243e3_dp*t47*t167 & 1394 *t331 - 0.80e2_dp/0.729e3_dp*t57*t178*t335 1395 t340 = t23*t338*t83 1396 t356 = -t67*t5*t319/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t70*t149*t323 & 1397 - 0.4e1_dp/0.27e2_dp*t73*t39*t327 - 0.32e2_dp/0.243e3_dp*t76*t167 & 1398 *t331 - 0.80e2_dp/0.729e3_dp*t79*t178*t335 1399 t358 = t66*t193*t356 1400 t359 = t3*t5 1401 t361 = t270*t103 1402 t363 = t359*t319*t361/0.9e1_dp 1403 t364 = t340*t103 1404 t366 = t102*R*t356 1405 t367 = t225*t366 1406 t388 = t232*t93*t97*t132*t3*t33*t314*t141/0.3e1_dp + 0.3e1_dp & 1407 *t240*t97*t10*t338 - 0.3e1_dp*t247*t97*t10*t356 + & 1408 t94*(-t317 - t340 + t358)*t95*t235 1409 t390 = t231*R*t388 1410 t392 = t84*t390/0.3e1_dp 1411 t394 = (t317 + t340 - t358 - t363 - t364 + t367 + t392)*t106 1412 t397 = t106*br_a1 1413 t399 = t2*t5*t33 1414 t403 = t338*t83 1415 t405 = t274*t356 1416 t409 = t397*t2 1417 t410 = t312*gamma 1418 t414 = ndrho*t141*t65*t114 1419 t423 = (-t317 - t340 + t358 - t363 - t364 + t367 + t392)*t117 1420 t426 = t117*br_a1 1421 t434 = t426*t2 1422 t443 = -0.2e1_dp*t394 + t394*t84 + t397*t399*t315/0.9e1_dp + t108 & 1423 *t403 - t108*t405 - t394*t66*t114 - t409*t410*t414/ & 1424 0.9e1_dp - t108*t338*t114 + t283*t366 + t110*t390/0.3e1_dp + & 1425 0.2e1_dp*t423 + t423*t84 + t426*t399*t315/0.9e1_dp + t119*t403 & 1426 - t119*t405 + t423*t66*t114 + t434*t410*t414/0.9e1_dp & 1427 + t119*t338*t114 - t301*t366 - t120*t390/0.3e1_dp 1428 e_ndrho = e_ndrho + (rho*t443*t102/0.8e1_dp - t124*t231*t388/0.24e2_dp)*sx 1429 t450 = t6*t33 1430 t455 = 0.4e1_dp/0.9e1_dp*t3*t450*gamma*t141*t109 1431 t456 = t450*gamma 1432 t459 = t31*t42 1433 t460 = t459*gamma 1434 t463 = t40*t51 1435 t464 = t463*gamma 1436 t467 = t49*t61 1437 t468 = t467*gamma 1438 t471 = t59*t185 1439 t472 = t471*gamma 1440 t475 = 0.4e1_dp/0.9e1_dp*t24*t456 + 0.16e2_dp/0.27e2_dp*t28*t460 + & 1441 0.16e2_dp/0.27e2_dp*t38*t464 + 0.128e3_dp/0.243e3_dp*t47*t468 + 0.320e3_dp & 1442 /0.729e3_dp*t57*t472 1443 t477 = t23*t475*t83 1444 t488 = 0.4e1_dp/0.9e1_dp*t67*t456 + 0.16e2_dp/0.27e2_dp*t70*t460 + & 1445 0.16e2_dp/0.27e2_dp*t73*t464 + 0.128e3_dp/0.243e3_dp*t76*t468 + 0.320e3_dp & 1446 /0.729e3_dp*t79*t472 1447 t490 = t66*t193*t488 1448 t493 = 0.4e1_dp/0.9e1_dp*t3*t456*t361 1449 t494 = t477*t103 1450 t496 = t102*R*t488 1451 t497 = t225*t496 1452 t499 = t232*t233*t96 1453 t516 = -0.4e1_dp/0.3e1_dp*t499*t359*t134*t141 + 0.3e1_dp*t240* & 1454 t97*t10*t475 - 0.3e1_dp*t247*t97*t10*t488 + t94*(t455 - & 1455 t477 + t490)*t95*t235 1456 t518 = t231*R*t516 1457 t520 = t84*t518/0.3e1_dp 1458 t522 = (-t455 + t477 - t490 + t493 - t494 + t497 + t520)*t106 1459 t525 = t2*t6 1460 t526 = t397*t525 1461 t527 = t134*t270 1462 t530 = t475*t83 1463 t532 = t274*t488 1464 t545 = (t455 - t477 + t490 + t493 - t494 + t497 + t520)*t117 1465 t548 = t426*t525 1466 t563 = -0.2e1_dp*t522 + t522*t84 - 0.4e1_dp/0.9e1_dp*t526*t527 + t108 & 1467 *t530 - t108*t532 - t522*t66*t114 + 0.4e1_dp/0.9e1_dp*t409 & 1468 *t456*t361 - t108*t475*t114 + t283*t496 + t110*t518/ & 1469 0.3e1_dp + 0.2e1_dp*t545 + t545*t84 - 0.4e1_dp/0.9e1_dp*t548*t527 + & 1470 t119*t530 - t119*t532 + t545*t66*t114 - 0.4e1_dp/0.9e1_dp*t434 & 1471 *t456*t361 + t119*t475*t114 - t301*t496 - t120*t518 & 1472 /0.3e1_dp 1473 e_tau = e_tau + (rho*t563*t102/0.8e1_dp - t124*t231*t516/0.24e2_dp)*sx 1474 t572 = t33*t141*t109 1475 t574 = t3*t6*t572/0.9e1_dp 1476 t585 = -t24*t450/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t28*t459 - 0.4e1_dp/ & 1477 0.27e2_dp*t38*t463 - 0.32e2_dp/0.243e3_dp*t47*t467 - 0.80e2_dp/0.729e3_dp & 1478 *t57*t471 1479 t587 = t23*t585*t83 1480 t598 = -t67*t450/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t70*t459 - 0.4e1_dp/ & 1481 0.27e2_dp*t73*t463 - 0.32e2_dp/0.243e3_dp*t76*t467 - 0.80e2_dp/0.729e3_dp & 1482 *t79*t471 1483 t600 = t66*t193*t598 1484 t605 = t3*t450*t141*t109*t103/0.9e1_dp 1485 t606 = t587*t103 1486 t608 = t102*R*t598 1487 t609 = t225*t608 1488 t628 = t499*t5*br_a1*t2*t33*t141/0.3e1_dp + 0.3e1_dp*t240* & 1489 t97*t10*t585 - 0.3e1_dp*t247*t97*t10*t598 + t94*(-t574 & 1490 - t587 + t600)*t95*t235 1491 t630 = t231*R*t628 1492 t632 = t84*t630/0.3e1_dp 1493 t634 = (t574 + t587 - t600 - t605 - t606 + t609 + t632)*t106 1494 t639 = t585*t83 1495 t641 = t274*t598 1496 t645 = t525*t33 1497 t655 = (-t574 - t587 + t600 - t605 - t606 + t609 + t632)*t117 1498 t672 = -0.2e1_dp*t634 + t634*t84 + t526*t572/0.9e1_dp + t108*t639 & 1499 - t108*t641 - t634*t66*t114 - t397*t645*t361/0.9e1_dp & 1500 - t108*t585*t114 + t283*t608 + t110*t630/0.3e1_dp + 0.2e1_dp* & 1501 t655 + t655*t84 + t548*t572/0.9e1_dp + t119*t639 - t119*t641 & 1502 + t655*t66*t114 + t426*t645*t361/0.9e1_dp + t119*t585 & 1503 *t114 - t301*t608 - t120*t630/0.3e1_dp 1504 e_laplace_rho = e_laplace_rho + (rho*t672*t102/0.8e1_dp - t124*t231*t628/0.24e2_dp)*sx 1505 END IF 1506 1507 END SUBROUTINE x_br_lsd_y_lte_0_cutoff_R_gt_b 1508 1509! ************************************************************************************************** 1510!> \brief Truncated long range part for y <= 0 and R <= b 1511!> \param rho ... 1512!> \param ndrho ... 1513!> \param tau ... 1514!> \param laplace_rho ... 1515!> \param e_0 ... 1516!> \param e_rho ... 1517!> \param e_ndrho ... 1518!> \param e_tau ... 1519!> \param e_laplace_rho ... 1520!> \param sx ... 1521!> \param R ... 1522!> \param gamma ... 1523!> \param grad_deriv ... 1524!> \par History 1525!> 12.2008 created [mguidon] 1526!> \author mguidon 1527! ************************************************************************************************** 1528 SUBROUTINE x_br_lsd_y_lte_0_cutoff_R_lte_b(rho, ndrho, tau, laplace_rho, e_0, & 1529 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1530 sx, R, gamma, grad_deriv) 1531 REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho 1532 REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho 1533 REAL(dp), INTENT(IN) :: sx, R, gamma 1534 INTEGER, INTENT(IN) :: grad_deriv 1535 1536 CHARACTER(len=*), PARAMETER :: routineN = 'x_br_lsd_y_lte_0_cutoff_R_lte_b', & 1537 routineP = moduleN//':'//routineN 1538 1539 REAL(KIND=dp) :: t1, t10, t100, t101, t102, t103, t104, t106, t108, t109, t110, t111, t113, & 1540 t115, t116, t118, t119, t121, t123, t128, t131, t133, t134, t137, t138, t140, t141, t143, & 1541 t146, t153, t154, t157, t159, t16, t160, t163, t167, t168, t17, t171, t172, t175, t177, & 1542 t178, t18, t182, t183, t186, t189, t190, t191, t194, t196, t197, t199, t2, t200, t21, & 1543 t22, t226, t229, t23, t232, t233, t234, t235, t237, t24, t242, t247, t249, t254, t256, & 1544 t264, t267, t27, t270, t272, t277, t28, t280, t281, t29, t292, t295, t298, t299, t3, t30, & 1545 t302, t31, t310, t314, t315, t317, t318, t32, t324, t328, t33 1546 REAL(KIND=dp) :: t332, t336, t339, t34, t341, t342, t359, t362, t368, t369, t37, t38, t383, & 1547 t385, t387, t39, t392, t395, t398, t4, t40, t402, t404, t409, t42, t420, t427, t428, & 1548 t429, t43, t432, t443, t444, t446, t450, t451, t454, t455, t458, t459, t46, t462, t463, & 1549 t466, t468, t469, t47, t48, t481, t484, t487, t49, t5, t50, t501, t504, t506, t51, t511, & 1550 t514, t517, t52, t521, t525, t529, t540, t547, t548, t549, t552, t56, t566, t57, t578, & 1551 t58, t580, t581, t59, t593, t596, t6, t61, t612, t613, t614, t616, t618, t62, t623, t626, & 1552 t629, t638, t65, t654, t655, t656, t659, t66, t67, t70, t73, t76 1553 REAL(KIND=dp) :: t79, t82, t83, t84, t85, t86, t87, t88, t89, t9, t90, t91, t93, t94, t95, & 1554 t96, t97, t99 1555 1556 IF (grad_deriv >= 0) THEN 1557 t1 = pi**(0.1e1_dp/0.3e1_dp) 1558 t2 = t1**2 1559 t3 = br_a1*t2 1560 t4 = rho**(0.1e1_dp/0.3e1_dp) 1561 t5 = t4**2 1562 t6 = t5*rho 1563 t9 = ndrho**2 1564 t10 = 0.1e1_dp/rho 1565 t16 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t9*t10/0.4e1_dp)/0.3e1_dp 1566 t17 = 0.1e1_dp/t16 1567 t18 = t6*t17 1568 t21 = 0.2e1_dp/0.3e1_dp*t3*t18 + br_a2 1569 t22 = ATAN(t21) 1570 t23 = -t22 + br_a3 1571 t24 = br_c1*t2 1572 t27 = t1*pi 1573 t28 = br_c2*t27 1574 t29 = rho**2 1575 t30 = t29*rho 1576 t31 = t4*t30 1577 t32 = t16**2 1578 t33 = 0.1e1_dp/t32 1579 t34 = t31*t33 1580 t37 = pi**2 1581 t38 = br_c3*t37 1582 t39 = t29**2 1583 t40 = t39*rho 1584 t42 = 0.1e1_dp/t32/t16 1585 t43 = t40*t42 1586 t46 = t2*t37 1587 t47 = br_c4*t46 1588 t48 = t39*t29 1589 t49 = t5*t48 1590 t50 = t32**2 1591 t51 = 0.1e1_dp/t50 1592 t52 = t49*t51 1593 t56 = t1*t37*pi 1594 t57 = br_c5*t56 1595 t58 = t39**2 1596 t59 = t4*t58 1597 t61 = 0.1e1_dp/t50/t16 1598 t62 = t59*t61 1599 t65 = br_c0 + 0.2e1_dp/0.3e1_dp*t24*t18 + 0.4e1_dp/0.9e1_dp*t28*t34 & 1600 + 0.8e1_dp/0.27e2_dp*t38*t43 + 0.16e2_dp/0.81e2_dp*t47*t52 + 0.32e2_dp & 1601 /0.243e3_dp*t57*t62 1602 t66 = t23*t65 1603 t67 = br_b1*t2 1604 t70 = br_b2*t27 1605 t73 = br_b3*t37 1606 t76 = br_b4*t46 1607 t79 = br_b5*t56 1608 t82 = br_b0 + 0.2e1_dp/0.3e1_dp*t67*t18 + 0.4e1_dp/0.9e1_dp*t70*t34 & 1609 + 0.8e1_dp/0.27e2_dp*t73*t43 + 0.16e2_dp/0.81e2_dp*t76*t52 + 0.32e2_dp & 1610 /0.243e3_dp*t79*t62 1611 t83 = 0.1e1_dp/t82 1612 t84 = t66*t83 1613 t85 = 8**(0.1e1_dp/0.3e1_dp) 1614 t86 = t23**2 1615 t87 = t86*t23 1616 t88 = t65**2 1617 t89 = t88*t65 1618 t90 = t87*t89 1619 t91 = t82**2 1620 t93 = 0.1e1_dp/t91/t82 1621 t94 = t90*t93 1622 t95 = EXP(-t84) 1623 t96 = 0.1e1_dp/0.3141592654e1_dp 1624 t97 = t95*t96 1625 t99 = t94*t97*t10 1626 t100 = t99**(0.1e1_dp/0.3e1_dp) 1627 t101 = 0.1e1_dp/t100 1628 t102 = REAL(t85, KIND=dp)*t101 1629 t103 = t102*R 1630 t104 = t84*t103 1631 t106 = EXP(0.2e1_dp*t104) 1632 t108 = t106*R 1633 t109 = t108*t23 1634 t110 = t65*t83 1635 t111 = t110*t102 1636 t113 = t106*t23 1637 t115 = t84 + t104 1638 t116 = EXP(t115) 1639 t118 = 0.2e1_dp*t106 - t109*t111 + t113*t110 + 0.2e1_dp + t84 + t104 & 1640 - 0.4e1_dp*t116 1641 t119 = rho*t118 1642 t121 = EXP(-t115) 1643 t123 = t121*REAL(t85, KIND=dp)*t101 1644 e_0 = e_0 + (t119*t123/0.8e1_dp)*sx 1645 END IF 1646 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 1647 t128 = t5*t17 1648 t131 = 0.1e1_dp/t4 1649 t133 = t33*gamma 1650 t134 = t133*t9 1651 t137 = 0.10e2_dp/0.9e1_dp*t3*t128 + t3*t131*t134/0.18e2_dp 1652 t138 = t21**2 1653 t140 = 0.1e1_dp/(0.1e1_dp + t138) 1654 t141 = t137*t140 1655 t143 = t83*REAL(t85, KIND=dp) 1656 t146 = t141*t65*t143*t101*R 1657 t153 = t4*t29 1658 t154 = t153*t33 1659 t157 = t4*rho 1660 t159 = t42*gamma 1661 t160 = t159*t9 1662 t163 = t39*t42 1663 t167 = t51*gamma 1664 t168 = t167*t9 1665 t171 = t5*t40 1666 t172 = t171*t51 1667 t175 = t5*t39 1668 t177 = t61*gamma 1669 t178 = t177*t9 1670 t182 = t4*t39*t30 1671 t183 = t182*t61 1672 t186 = t4*t48 1673 t189 = 0.1e1_dp/t50/t32 1674 t190 = t189*gamma 1675 t191 = t190*t9 1676 t194 = 0.10e2_dp/0.9e1_dp*t24*t128 + t24*t131*t134/0.18e2_dp + 0.40e2_dp & 1677 /0.27e2_dp*t28*t154 + 0.2e1_dp/0.27e2_dp*t28*t157*t160 + & 1678 0.40e2_dp/0.27e2_dp*t38*t163 + 0.2e1_dp/0.27e2_dp*t38*t30*t168 & 1679 + 0.320e3_dp/0.243e3_dp*t47*t172 + 0.16e2_dp/0.243e3_dp*t47*t175* & 1680 t178 + 0.800e3_dp/0.729e3_dp*t57*t183 + 0.40e2_dp/0.729e3_dp*t57* & 1681 t186*t191 1682 t196 = t23*t194*t83 1683 t197 = t196*t103 1684 t199 = 0.1e1_dp/t91 1685 t200 = t66*t199 1686 t226 = 0.10e2_dp/0.9e1_dp*t67*t128 + t67*t131*t134/0.18e2_dp + 0.40e2_dp & 1687 /0.27e2_dp*t70*t154 + 0.2e1_dp/0.27e2_dp*t70*t157*t160 + & 1688 0.40e2_dp/0.27e2_dp*t73*t163 + 0.2e1_dp/0.27e2_dp*t73*t30*t168 & 1689 + 0.320e3_dp/0.243e3_dp*t76*t172 + 0.16e2_dp/0.243e3_dp*t76*t175* & 1690 t178 + 0.800e3_dp/0.729e3_dp*t79*t183 + 0.40e2_dp/0.729e3_dp*t79* & 1691 t186*t191 1692 t229 = t200*t102*R*t226 1693 t232 = 0.1e1_dp/t100/t99 1694 t233 = REAL(t85, KIND=dp)*t232 1695 t234 = t86*t89 1696 t235 = t93*t95 1697 t237 = t96*t10 1698 t242 = t87*t88*t93 1699 t247 = t91**2 1700 t249 = t90/t247 1701 t254 = t141*t110 1702 t256 = t66*t199*t226 1703 t264 = -0.3e1_dp*t234*t235*t237*t141 + 0.3e1_dp*t242*t97*t10 & 1704 *t194 - 0.3e1_dp*t249*t97*t10*t226 + t94*(t254 - t196 + & 1705 t256)*t95*t237 - t94*t97/t29 1706 t267 = t84*t233*R*t264 1707 t270 = (-0.2e1_dp*t146 + 0.2e1_dp*t197 - 0.2e1_dp*t229 - 0.2e1_dp/0.3e1_dp & 1708 *t267)*t106 1709 t272 = R*t23 1710 t277 = t194*t83 1711 t280 = t108*t66 1712 t281 = t199*REAL(t85, KIND=dp) 1713 t292 = t140*t65*t83 1714 t295 = t65*t199 1715 t298 = t267/0.3e1_dp 1716 t299 = -t254 + t196 - t256 - t146 + t197 - t229 - t298 1717 t302 = 0.2e1_dp*t270 - t270*t272*t111 + t108*t141*t111 - t109 & 1718 *t277*t102 + t280*t281*t101*t226 + t280*t143*t232* & 1719 t264/0.3e1_dp + t270*t84 - t106*t137*t292 + t113*t277 - t113 & 1720 *t295*t226 - t254 + t196 - t256 - t146 + t197 - t229 - t298 & 1721 - 0.4e1_dp*t299*t116 1722 t310 = t119*t121 1723 e_rho = e_rho + (t118*t121*t102/0.8e1_dp + rho*t302*t123/0.8e1_dp - t119 & 1724 *t299*t123/0.8e1_dp - t310*t233*t264/0.24e2_dp)*sx 1725 t314 = t3*t5 1726 t315 = t133*ndrho 1727 t317 = t292*t103 1728 t318 = t314*t315*t317 1729 t324 = t159*ndrho 1730 t328 = t167*ndrho 1731 t332 = t177*ndrho 1732 t336 = t190*ndrho 1733 t339 = -t24*t5*t315/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t28*t153*t324 & 1734 - 0.4e1_dp/0.27e2_dp*t38*t39*t328 - 0.32e2_dp/0.243e3_dp*t47*t171 & 1735 *t332 - 0.80e2_dp/0.729e3_dp*t57*t182*t336 1736 t341 = t23*t339*t83 1737 t342 = t341*t103 1738 t359 = -t67*t5*t315/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t70*t153*t324 & 1739 - 0.4e1_dp/0.27e2_dp*t73*t39*t328 - 0.32e2_dp/0.243e3_dp*t76*t171 & 1740 *t332 - 0.80e2_dp/0.729e3_dp*t79*t182*t336 1741 t362 = t200*t102*R*t359 1742 t368 = gamma*ndrho 1743 t369 = t368*t140 1744 t383 = t368*t292 1745 t385 = t3*t5*t33*t383/0.9e1_dp 1746 t387 = t66*t199*t359 1747 t392 = t234*t93*t97*t131*t3*t33*t369/0.3e1_dp + 0.3e1_dp* & 1748 t242*t97*t10*t339 - 0.3e1_dp*t249*t97*t10*t359 + t94* & 1749 (-t385 - t341 + t387)*t95*t237 1750 t395 = t84*t233*R*t392 1751 t398 = (0.2e1_dp/0.9e1_dp*t318 + 0.2e1_dp*t342 - 0.2e1_dp*t362 - 0.2e1_dp & 1752 /0.3e1_dp*t395)*t106 1753 t402 = t108*br_a1 1754 t404 = t2*t5*t33 1755 t409 = t339*t83 1756 t420 = t106*br_a1 1757 t427 = t318/0.9e1_dp 1758 t428 = t395/0.3e1_dp 1759 t429 = t385 + t341 - t387 + t427 + t342 - t362 - t428 1760 t432 = 0.2e1_dp*t398 - t398*t272*t111 - t402*t404*t369*t111 & 1761 /0.9e1_dp - t109*t409*t102 + t280*t281*t101*t359 + t280 & 1762 *t143*t232*t392/0.3e1_dp + t398*t84 + t420*t404*t383/0.9e1_dp & 1763 + t113*t409 - t113*t295*t359 + t385 + t341 - t387 + t427 & 1764 + t342 - t362 - t428 - 0.4e1_dp*t429*t116 1765 e_ndrho = e_ndrho + (rho*t432*t123/0.8e1_dp - t119*t429*t123/0.8e1_dp - t310 & 1766 *t233*t392/0.24e2_dp)*sx 1767 t443 = t6*t33 1768 t444 = t443*gamma 1769 t446 = t3*t444*t317 1770 t450 = t31*t42 1771 t451 = t450*gamma 1772 t454 = t40*t51 1773 t455 = t454*gamma 1774 t458 = t49*t61 1775 t459 = t458*gamma 1776 t462 = t59*t189 1777 t463 = t462*gamma 1778 t466 = 0.4e1_dp/0.9e1_dp*t24*t444 + 0.16e2_dp/0.27e2_dp*t28*t451 + & 1779 0.16e2_dp/0.27e2_dp*t38*t455 + 0.128e3_dp/0.243e3_dp*t47*t459 + 0.320e3_dp & 1780 /0.729e3_dp*t57*t463 1781 t468 = t23*t466*t83 1782 t469 = t468*t103 1783 t481 = 0.4e1_dp/0.9e1_dp*t67*t444 + 0.16e2_dp/0.27e2_dp*t70*t451 + & 1784 0.16e2_dp/0.27e2_dp*t73*t455 + 0.128e3_dp/0.243e3_dp*t76*t459 + 0.320e3_dp & 1785 /0.729e3_dp*t79*t463 1786 t484 = t200*t102*R*t481 1787 t487 = t234*t235*t96 1788 t501 = gamma*t140 1789 t504 = 0.4e1_dp/0.9e1_dp*t3*t443*t501*t110 1790 t506 = t66*t199*t481 1791 t511 = -0.4e1_dp/0.3e1_dp*t487*t314*t133*t140 + 0.3e1_dp*t242* & 1792 t97*t10*t466 - 0.3e1_dp*t249*t97*t10*t481 + t94*(t504 - & 1793 t468 + t506)*t95*t237 1794 t514 = t84*t233*R*t511 1795 t517 = (-0.8e1_dp/0.9e1_dp*t446 + 0.2e1_dp*t469 - 0.2e1_dp*t484 - 0.2e1_dp & 1796 /0.3e1_dp*t514)*t106 1797 t521 = t2*t6 1798 t525 = t143*t101 1799 t529 = t466*t83 1800 t540 = t420*t521 1801 t547 = 0.4e1_dp/0.9e1_dp*t446 1802 t548 = t514/0.3e1_dp 1803 t549 = -t504 + t468 - t506 - t547 + t469 - t484 - t548 1804 t552 = 0.2e1_dp*t517 - t517*t272*t111 + 0.4e1_dp/0.9e1_dp*t402*t521 & 1805 *t33*t501*t65*t525 - t109*t529*t102 + t280*t281* & 1806 t101*t481 + t280*t143*t232*t511/0.3e1_dp + t517*t84 - 0.4e1_dp & 1807 /0.9e1_dp*t540*t133*t292 + t113*t529 - t113*t295*t481 & 1808 - t504 + t468 - t506 - t547 + t469 - t484 - t548 - 0.4e1_dp*t549 & 1809 *t116 1810 e_tau = e_tau + (rho*t552*t123/0.8e1_dp - t119*t549*t123/0.8e1_dp - t310 & 1811 *t233*t511/0.24e2_dp)*sx 1812 t566 = t3*t443*t140*t110*t103 1813 t578 = -t24*t443/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t28*t450 - 0.4e1_dp/ & 1814 0.27e2_dp*t38*t454 - 0.32e2_dp/0.243e3_dp*t47*t458 - 0.80e2_dp/0.729e3_dp & 1815 *t57*t462 1816 t580 = t23*t578*t83 1817 t581 = t580*t103 1818 t593 = -t67*t443/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t70*t450 - 0.4e1_dp/ & 1819 0.27e2_dp*t73*t454 - 0.32e2_dp/0.243e3_dp*t76*t458 - 0.80e2_dp/0.729e3_dp & 1820 *t79*t462 1821 t596 = t200*t102*R*t593 1822 t612 = t3*t6 1823 t613 = t33*t140 1824 t614 = t613*t110 1825 t616 = t612*t614/0.9e1_dp 1826 t618 = t66*t199*t593 1827 t623 = t487*t5*br_a1*t2*t33*t140/0.3e1_dp + 0.3e1_dp*t242* & 1828 t97*t10*t578 - 0.3e1_dp*t249*t97*t10*t593 + t94*(-t616 & 1829 - t580 + t618)*t95*t237 1830 t626 = t84*t233*R*t623 1831 t629 = (0.2e1_dp/0.9e1_dp*t566 + 0.2e1_dp*t581 - 0.2e1_dp*t596 - 0.2e1_dp & 1832 /0.3e1_dp*t626)*t106 1833 t638 = t578*t83 1834 t654 = t566/0.9e1_dp 1835 t655 = t626/0.3e1_dp 1836 t656 = t616 + t580 - t618 + t654 + t581 - t596 - t655 1837 t659 = 0.2e1_dp*t629 - t629*t272*t111 - t108*t612*t613*t65 & 1838 *t525/0.9e1_dp - t109*t638*t102 + t280*t281*t101*t593 + & 1839 t280*t143*t232*t623/0.3e1_dp + t629*t84 + t540*t614/0.9e1_dp & 1840 + t113*t638 - t113*t295*t593 + t616 + t580 - t618 + t654 & 1841 + t581 - t596 - t655 - 0.4e1_dp*t656*t116 1842 e_laplace_rho = e_laplace_rho + (rho*t659*t123/0.8e1_dp - t119*t656*t123/0.8e1_dp - t310 & 1843 *t233*t623/0.24e2_dp)*sx 1844 END IF 1845 1846 END SUBROUTINE x_br_lsd_y_lte_0_cutoff_R_lte_b 1847 1848! ************************************************************************************************** 1849!> \brief Truncated long range part for y > 0 1850!> \param rho ... 1851!> \param ndrho ... 1852!> \param tau ... 1853!> \param laplace_rho ... 1854!> \param e_0 ... 1855!> \param e_rho ... 1856!> \param e_ndrho ... 1857!> \param e_tau ... 1858!> \param e_laplace_rho ... 1859!> \param sx ... 1860!> \param R ... 1861!> \param gamma ... 1862!> \param grad_deriv ... 1863!> \par History 1864!> 12.2008 created [mguidon] 1865!> \author mguidon 1866! ************************************************************************************************** 1867 SUBROUTINE x_br_lsd_y_gt_0_cutoff(rho, ndrho, tau, laplace_rho, e_0, & 1868 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1869 sx, R, gamma, grad_deriv) 1870 REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho 1871 REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho 1872 REAL(dp), INTENT(IN) :: sx, R, gamma 1873 INTEGER, INTENT(IN) :: grad_deriv 1874 1875 CHARACTER(len=*), PARAMETER :: routineN = 'x_br_lsd_y_gt_0_cutoff', & 1876 routineP = moduleN//':'//routineN 1877 1878 REAL(KIND=dp) :: brval, t1, t10, t103, t104, t11, t111, t116, t14, t15, t2, t21, t25, t26, & 1879 t28, t3, t31, t33, t37, t4, t44, t45, t46, t5, t50, t56, t58, t6, t62, t65, t69, t71, & 1880 t75, t77, t8, t81, t84, t85, t9 1881 1882 t1 = 8**(0.1e1_dp/0.3e1_dp) 1883 t2 = t1**2 1884 t3 = 0.1e1_dp/br_BB 1885 t4 = pi**(0.1e1_dp/0.3e1_dp) 1886 t5 = t4**2 1887 t6 = 0.1e1_dp/t5 1888 t8 = rho**(0.1e1_dp/0.3e1_dp) 1889 t9 = t8**2 1890 t10 = t9*rho 1891 t11 = 0.1e1_dp/t10 1892 t14 = ndrho**2 1893 t15 = 0.1e1_dp/rho 1894 t21 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t14*t15/0.4e1_dp)/0.3e1_dp 1895 t25 = br_BB**2 1896 t26 = t4*pi 1897 t28 = rho**2 1898 t31 = t21**2 1899 t33 = t8*t28*rho/t31 1900 t37 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t26*t33) 1901 t44 = LOG(0.1500000000000000e1_dp*t3*t6*t11*t21 + 0.3e1_dp/0.2e1_dp & 1902 *t37*t3*t6*t11*t21) 1903 t45 = t44 + 0.2e1_dp 1904 t46 = t45**2 1905 t50 = t10/t21 1906 t56 = pi**2 1907 t58 = t28**2 1908 t62 = t58*rho/t31/t21 1909 t65 = t5*t56 1910 t69 = t31**2 1911 t71 = t9*t58*t28/t69 1912 t75 = t4*t56*pi 1913 t77 = t58**2 1914 t81 = t8*t77/t69/t21 1915 t84 = br_d0 + 0.2e1_dp/0.3e1_dp*br_d1*t5*t50 + 0.4e1_dp/0.9e1_dp*br_d2 & 1916 *t26*t33 + 0.8e1_dp/0.27e2_dp*br_d3*t56*t62 + 0.16e2_dp/0.81e2_dp & 1917 *br_d4*t65*t71 + 0.32e2_dp/0.243e3_dp*br_d5*t75*t81 1918 t85 = t84**2 1919 t103 = br_e0 + 0.2e1_dp/0.3e1_dp*br_e1*t5*t50 + 0.4e1_dp/0.9e1_dp*br_e2 & 1920 *t26*t33 + 0.8e1_dp/0.27e2_dp*br_e3*t56*t62 + 0.16e2_dp/0.81e2_dp & 1921 *br_e4*t65*t71 + 0.32e2_dp/0.243e3_dp*br_e5*t75*t81 1922 t104 = t103**2 1923 t111 = EXP(-t45*t84/t103) 1924 t116 = (t46*t45*t85*t84/t104/t103*t111/0.3141592654e1_dp & 1925 *t15)**(0.1e1_dp/0.3e1_dp) 1926 brval = REAL(t2, KIND=dp)*t116/0.8e1_dp 1927 1928 IF (R > brval) THEN 1929 CALL x_br_lsd_y_gt_0_cutoff_R_gt_b(rho, ndrho, tau, laplace_rho, e_0, & 1930 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1931 sx, R, gamma, grad_deriv) 1932 ELSE 1933 CALL x_br_lsd_y_gt_0_cutoff_R_lte_b(rho, ndrho, tau, laplace_rho, e_0, & 1934 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1935 sx, R, gamma, grad_deriv) 1936 END IF 1937 1938 END SUBROUTINE x_br_lsd_y_gt_0_cutoff 1939 1940! ************************************************************************************************** 1941!> \brief Truncated long range part for y > 0 and R > b 1942!> \param rho ... 1943!> \param ndrho ... 1944!> \param tau ... 1945!> \param laplace_rho ... 1946!> \param e_0 ... 1947!> \param e_rho ... 1948!> \param e_ndrho ... 1949!> \param e_tau ... 1950!> \param e_laplace_rho ... 1951!> \param sx ... 1952!> \param R ... 1953!> \param gamma ... 1954!> \param grad_deriv ... 1955!> \par History 1956!> 12.2008 created [mguidon] 1957!> \author mguidon 1958! ************************************************************************************************** 1959 SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_gt_b(rho, ndrho, tau, laplace_rho, e_0, & 1960 e_rho, e_ndrho, e_tau, e_laplace_rho, & 1961 sx, R, gamma, grad_deriv) 1962 1963 REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho 1964 REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho 1965 REAL(dp), INTENT(IN) :: sx, R, gamma 1966 INTEGER, INTENT(IN) :: grad_deriv 1967 1968 CHARACTER(len=*), PARAMETER :: routineN = 'x_br_lsd_y_gt_0_cutoff_R_gt_b', & 1969 routineP = moduleN//':'//routineN 1970 1971 REAL(KIND=dp) :: t1, t100, t101, t102, t103, t104, t105, t106, t108, t109, t110, t111, t112, & 1972 t114, t115, t116, t117, t118, t119, t12, t121, t123, t124, t125, t129, t13, t132, t134, & 1973 t135, t138, t139, t145, t152, t155, t158, t159, t162, t164, t165, t176, t179, t180, t181, & 1974 t182, t183, t186, t188, t189, t19, t197, t2, t20, t201, t202, t205, t206, t209, t211, & 1975 t212, t216, t217, t220, t223, t224, t225, t228, t23, t230, t231, t24, t25, t257, t259, & 1976 t26, t261, t262, t263, t265, t266, t269, t27, t272, t273, t278, t28, t283, t285, t29, & 1977 t297, t299, t3, t30, t301, t303, t306, t307, t308, t31, t310, t312 1978 REAL(KIND=dp) :: t313, t321, t326, t329, t339, t343, t35, t351, t354, t355, t36, t363, t364, & 1979 t365, t367, t37, t371, t375, t379, t383, t386, t388, t4, t404, t406, t408, t409, t41, & 1980 t411, t412, t42, t428, t43, t430, t432, t434, t437, t439, t44, t441, t45, t453, t456, & 1981 t46, t469, t479, t480, t485, t486, t487, t49, t490, t491, t494, t495, t498, t499, t5, & 1982 t502, t503, t506, t508, t519, t52, t521, t523, t524, t526, t527, t53, t54, t543, t545, & 1983 t547, t549, t55, t552, t554, t556, t568, t57, t571, t58, t584, t599, t6, t600, t601, t61, & 1984 t612, t614, t62, t625, t627, t629, t63, t630, t632, t633, t64, t649 1985 REAL(KIND=dp) :: t65, t651, t653, t655, t658, t66, t660, t662, t67, t674, t677, t690, t7, & 1986 t71, t72, t73, t74, t76, t77, t8, t80, t81, t82, t85, t88, t9, t91, t94, t97, t98, t99 1987 1988 IF (grad_deriv >= 0) THEN 1989 t1 = 0.1e1_dp/br_BB 1990 t2 = pi**(0.1e1_dp/0.3e1_dp) 1991 t3 = t2**2 1992 t4 = 0.1e1_dp/t3 1993 t5 = t1*t4 1994 t6 = rho**(0.1e1_dp/0.3e1_dp) 1995 t7 = t6**2 1996 t8 = t7*rho 1997 t9 = 0.1e1_dp/t8 1998 t12 = ndrho**2 1999 t13 = 0.1e1_dp/rho 2000 t19 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t12*t13/0.4e1_dp)/0.3e1_dp 2001 t20 = t9*t19 2002 t23 = br_BB**2 2003 t24 = t2*pi 2004 t25 = t23*t24 2005 t26 = rho**2 2006 t27 = t26*rho 2007 t28 = t6*t27 2008 t29 = t19**2 2009 t30 = 0.1e1_dp/t29 2010 t31 = t28*t30 2011 t35 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t31) 2012 t36 = t35*t1 2013 t37 = t4*t9 2014 t41 = 0.1500000000000000e1_dp*t5*t20 + 0.3e1_dp/0.2e1_dp*t36*t37* & 2015 t19 2016 t42 = LOG(t41) 2017 t43 = t42 + 0.2e1_dp 2018 t44 = br_d1*t3 2019 t45 = 0.1e1_dp/t19 2020 t46 = t8*t45 2021 t49 = br_d2*t24 2022 t52 = pi**2 2023 t53 = br_d3*t52 2024 t54 = t26**2 2025 t55 = t54*rho 2026 t57 = 0.1e1_dp/t29/t19 2027 t58 = t55*t57 2028 t61 = t3*t52 2029 t62 = br_d4*t61 2030 t63 = t54*t26 2031 t64 = t7*t63 2032 t65 = t29**2 2033 t66 = 0.1e1_dp/t65 2034 t67 = t64*t66 2035 t71 = t2*t52*pi 2036 t72 = br_d5*t71 2037 t73 = t54**2 2038 t74 = t6*t73 2039 t76 = 0.1e1_dp/t65/t19 2040 t77 = t74*t76 2041 t80 = br_d0 + 0.2e1_dp/0.3e1_dp*t44*t46 + 0.4e1_dp/0.9e1_dp*t49*t31 & 2042 + 0.8e1_dp/0.27e2_dp*t53*t58 + 0.16e2_dp/0.81e2_dp*t62*t67 + 0.32e2_dp & 2043 /0.243e3_dp*t72*t77 2044 t81 = t43*t80 2045 t82 = br_e1*t3 2046 t85 = br_e2*t24 2047 t88 = br_e3*t52 2048 t91 = br_e4*t61 2049 t94 = br_e5*t71 2050 t97 = br_e0 + 0.2e1_dp/0.3e1_dp*t82*t46 + 0.4e1_dp/0.9e1_dp*t85*t31 & 2051 + 0.8e1_dp/0.27e2_dp*t88*t58 + 0.16e2_dp/0.81e2_dp*t91*t67 + 0.32e2_dp & 2052 /0.243e3_dp*t94*t77 2053 t98 = 0.1e1_dp/t97 2054 t99 = t81*t98 2055 t100 = 8**(0.1e1_dp/0.3e1_dp) 2056 t101 = t43**2 2057 t102 = t101*t43 2058 t103 = t80**2 2059 t104 = t103*t80 2060 t105 = t102*t104 2061 t106 = t97**2 2062 t108 = 0.1e1_dp/t106/t97 2063 t109 = t105*t108 2064 t110 = EXP(-t99) 2065 t111 = 0.1e1_dp/0.3141592654e1_dp 2066 t112 = t110*t111 2067 t114 = t109*t112*t13 2068 t115 = t114**(0.1e1_dp/0.3e1_dp) 2069 t116 = 0.1e1_dp/t115 2070 t117 = REAL(t100, KIND=dp)*t116 2071 t118 = t117*R 2072 t119 = t99*t118 2073 t121 = EXP(t99 - t119) 2074 t123 = t121*t43 2075 t124 = t80*t98 2076 t125 = t123*t124 2077 t129 = t98*REAL(t100, KIND=dp)*t116*R 2078 t132 = EXP(-t99 - t119) 2079 t134 = t132*t43 2080 t135 = t134*t124 2081 t138 = -0.2e1_dp*t121 + t125 - t123*t80*t129 + 0.2e1_dp*t132 + t135 & 2082 + t134*t80*t129 2083 t139 = rho*t138 2084 e_0 = e_0 + (t139*t117/0.8e1_dp)*sx 2085 END IF 2086 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 2087 t145 = 0.1e1_dp/t7/t26 2088 t152 = 0.1e1_dp/t7/t27*gamma*t12 2089 t155 = 0.1e1_dp/t35 2090 t158 = t6*t26 2091 t159 = t158*t30 2092 t162 = t6*rho 2093 t164 = t57*gamma 2094 t165 = t164*t12 2095 t176 = t36*t4 2096 t179 = -0.2500000000e1_dp*t5*t145*t19 - 0.1250000000e0_dp*t5*t152 & 2097 + 0.3e1_dp/0.4e1_dp*t155*t1*t4*t20*(0.40e2_dp/0.27e2_dp*t25 & 2098 *t159 + 0.2e1_dp/0.27e2_dp*t25*t162*t165) - 0.5e1_dp/0.2e1_dp*t36 & 2099 *t4*t145*t19 - t176*t152/0.8e1_dp 2100 t180 = 0.1e1_dp/t41 2101 t181 = t179*t180 2102 t182 = t181*t124 2103 t183 = t7*t45 2104 t186 = 0.1e1_dp/t6 2105 t188 = t30*gamma 2106 t189 = t188*t12 2107 t197 = t54*t57 2108 t201 = t66*gamma 2109 t202 = t201*t12 2110 t205 = t7*t55 2111 t206 = t205*t66 2112 t209 = t7*t54 2113 t211 = t76*gamma 2114 t212 = t211*t12 2115 t216 = t6*t54*t27 2116 t217 = t216*t76 2117 t220 = t6*t63 2118 t223 = 0.1e1_dp/t65/t29 2119 t224 = t223*gamma 2120 t225 = t224*t12 2121 t228 = 0.10e2_dp/0.9e1_dp*t44*t183 + t44*t186*t189/0.18e2_dp + 0.40e2_dp & 2122 /0.27e2_dp*t49*t159 + 0.2e1_dp/0.27e2_dp*t49*t162*t165 + & 2123 0.40e2_dp/0.27e2_dp*t53*t197 + 0.2e1_dp/0.27e2_dp*t53*t27*t202 & 2124 + 0.320e3_dp/0.243e3_dp*t62*t206 + 0.16e2_dp/0.243e3_dp*t62*t209* & 2125 t212 + 0.800e3_dp/0.729e3_dp*t72*t217 + 0.40e2_dp/0.729e3_dp*t72* & 2126 t220*t225 2127 t230 = t43*t228*t98 2128 t231 = 0.1e1_dp/t106 2129 t257 = 0.10e2_dp/0.9e1_dp*t82*t183 + t82*t186*t189/0.18e2_dp + 0.40e2_dp & 2130 /0.27e2_dp*t85*t159 + 0.2e1_dp/0.27e2_dp*t85*t162*t165 + & 2131 0.40e2_dp/0.27e2_dp*t88*t197 + 0.2e1_dp/0.27e2_dp*t88*t27*t202 & 2132 + 0.320e3_dp/0.243e3_dp*t91*t206 + 0.16e2_dp/0.243e3_dp*t91*t209* & 2133 t212 + 0.800e3_dp/0.729e3_dp*t94*t217 + 0.40e2_dp/0.729e3_dp*t94* & 2134 t220*t225 2135 t259 = t81*t231*t257 2136 t261 = t181*t80*t129 2137 t262 = t230*t118 2138 t263 = t81*t231 2139 t265 = t117*R*t257 2140 t266 = t263*t265 2141 t269 = REAL(t100, KIND=dp)/t115/t114 2142 t272 = t101*t104*t108*t110 2143 t273 = t111*t13 2144 t278 = t102*t103*t108 2145 t283 = t106**2 2146 t285 = t105/t283 2147 t297 = 0.3e1_dp*t272*t273*t181 + 0.3e1_dp*t278*t112*t13*t228 & 2148 - 0.3e1_dp*t285*t112*t13*t257 + t109*(-t182 - t230 + t259) & 2149 *t110*t273 - t109*t112/t26 2150 t299 = t269*R*t297 2151 t301 = t99*t299/0.3e1_dp 2152 t303 = (t182 + t230 - t259 - t261 - t262 + t266 + t301)*t121 2153 t306 = t121*t179 2154 t307 = t180*t80 2155 t308 = t307*t98 2156 t310 = t228*t98 2157 t312 = t80*t231 2158 t313 = t312*t257 2159 t321 = t123*t312 2160 t326 = (-t182 - t230 + t259 - t261 - t262 + t266 + t301)*t132 2161 t329 = t132*t179 2162 t339 = t134*t312 2163 t343 = -0.2e1_dp*t303 + t303*t99 + t306*t308 + t123*t310 - t123 & 2164 *t313 - t303*t81*t129 - t306*t307*t129 - t123*t228* & 2165 t129 + t321*t265 + t125*t299/0.3e1_dp + 0.2e1_dp*t326 + t326*t99 & 2166 + t329*t308 + t134*t310 - t134*t313 + t326*t81*t129 + & 2167 t329*t307*t129 + t134*t228*t129 - t339*t265 - t135*t299 & 2168 /0.3e1_dp 2169 e_rho = e_rho + (t138*REAL(t100, KIND=dp)*t116/0.8e1_dp + rho*t343*t117/0.8e1_dp & 2170 - t139*t269*t297/0.24e2_dp)*sx 2171 t351 = t145*gamma*ndrho 2172 t354 = t155*br_BB 2173 t355 = t354*t3 2174 t363 = 0.2500000000000000e0_dp*t5*t351 - t355*t7*t30*gamma*ndrho & 2175 /0.9e1_dp + t176*t351/0.4e1_dp 2176 t364 = t363*t180 2177 t365 = t364*t124 2178 t367 = t188*ndrho 2179 t371 = t164*ndrho 2180 t375 = t201*ndrho 2181 t379 = t211*ndrho 2182 t383 = t224*ndrho 2183 t386 = -t44*t7*t367/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t49*t158*t371 & 2184 - 0.4e1_dp/0.27e2_dp*t53*t54*t375 - 0.32e2_dp/0.243e3_dp*t62*t205 & 2185 *t379 - 0.80e2_dp/0.729e3_dp*t72*t216*t383 2186 t388 = t43*t386*t98 2187 t404 = -t82*t7*t367/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t85*t158*t371 & 2188 - 0.4e1_dp/0.27e2_dp*t88*t54*t375 - 0.32e2_dp/0.243e3_dp*t91*t205 & 2189 *t379 - 0.80e2_dp/0.729e3_dp*t94*t216*t383 2190 t406 = t81*t231*t404 2191 t408 = t364*t80*t129 2192 t409 = t388*t118 2193 t411 = t117*R*t404 2194 t412 = t263*t411 2195 t428 = 0.3e1_dp*t272*t273*t364 + 0.3e1_dp*t278*t112*t13*t386 & 2196 - 0.3e1_dp*t285*t112*t13*t404 + t109*(-t365 - t388 + t406) & 2197 *t110*t273 2198 t430 = t269*R*t428 2199 t432 = t99*t430/0.3e1_dp 2200 t434 = (t365 + t388 - t406 - t408 - t409 + t412 + t432)*t121 2201 t437 = t121*t363 2202 t439 = t386*t98 2203 t441 = t312*t404 2204 t453 = (-t365 - t388 + t406 - t408 - t409 + t412 + t432)*t132 2205 t456 = t132*t363 2206 t469 = -0.2e1_dp*t434 + t434*t99 + t437*t308 + t123*t439 - t123 & 2207 *t441 - t434*t81*t129 - t437*t307*t129 - t123*t386* & 2208 t129 + t321*t411 + t125*t430/0.3e1_dp + 0.2e1_dp*t453 + t453*t99 & 2209 + t456*t308 + t134*t439 - t134*t441 + t453*t81*t129 + & 2210 t456*t307*t129 + t134*t386*t129 - t339*t411 - t135*t430 & 2211 /0.3e1_dp 2212 e_ndrho = e_ndrho + (rho*t469*t117/0.8e1_dp - t139*t269*t428/0.24e2_dp)*sx 2213 t479 = t8*t30 2214 t480 = t479*gamma 2215 t485 = -0.1000000000e1_dp*t5*t9*gamma + 0.4e1_dp/0.9e1_dp*t355*t480 & 2216 - t36*t37*gamma 2217 t486 = t485*t180 2218 t487 = t486*t124 2219 t490 = t28*t57 2220 t491 = t490*gamma 2221 t494 = t55*t66 2222 t495 = t494*gamma 2223 t498 = t64*t76 2224 t499 = t498*gamma 2225 t502 = t74*t223 2226 t503 = t502*gamma 2227 t506 = 0.4e1_dp/0.9e1_dp*t44*t480 + 0.16e2_dp/0.27e2_dp*t49*t491 + & 2228 0.16e2_dp/0.27e2_dp*t53*t495 + 0.128e3_dp/0.243e3_dp*t62*t499 + 0.320e3_dp & 2229 /0.729e3_dp*t72*t503 2230 t508 = t43*t506*t98 2231 t519 = 0.4e1_dp/0.9e1_dp*t82*t480 + 0.16e2_dp/0.27e2_dp*t85*t491 + & 2232 0.16e2_dp/0.27e2_dp*t88*t495 + 0.128e3_dp/0.243e3_dp*t91*t499 + 0.320e3_dp & 2233 /0.729e3_dp*t94*t503 2234 t521 = t81*t231*t519 2235 t523 = t486*t80*t129 2236 t524 = t508*t118 2237 t526 = t117*R*t519 2238 t527 = t263*t526 2239 t543 = 0.3e1_dp*t272*t273*t486 + 0.3e1_dp*t278*t112*t13*t506 & 2240 - 0.3e1_dp*t285*t112*t13*t519 + t109*(-t487 - t508 + t521) & 2241 *t110*t273 2242 t545 = t269*R*t543 2243 t547 = t99*t545/0.3e1_dp 2244 t549 = (t487 + t508 - t521 - t523 - t524 + t527 + t547)*t121 2245 t552 = t121*t485 2246 t554 = t506*t98 2247 t556 = t312*t519 2248 t568 = (-t487 - t508 + t521 - t523 - t524 + t527 + t547)*t132 2249 t571 = t132*t485 2250 t584 = -0.2e1_dp*t549 + t549*t99 + t552*t308 + t123*t554 - t123 & 2251 *t556 - t549*t81*t129 - t552*t307*t129 - t123*t506* & 2252 t129 + t321*t526 + t125*t545/0.3e1_dp + 0.2e1_dp*t568 + t568*t99 & 2253 + t571*t308 + t134*t554 - t134*t556 + t568*t81*t129 + & 2254 t571*t307*t129 + t134*t506*t129 - t339*t526 - t135*t545 & 2255 /0.3e1_dp 2256 e_tau = e_tau + (rho*t584*t117/0.8e1_dp - t139*t269*t543/0.24e2_dp)*sx 2257 t599 = 0.2500000000000000e0_dp*t5*t9 - t354*t3*t8*t30/0.9e1_dp & 2258 + t36*t37/0.4e1_dp 2259 t600 = t599*t180 2260 t601 = t600*t124 2261 t612 = -t44*t479/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t49*t490 - 0.4e1_dp/ & 2262 0.27e2_dp*t53*t494 - 0.32e2_dp/0.243e3_dp*t62*t498 - 0.80e2_dp/0.729e3_dp & 2263 *t72*t502 2264 t614 = t43*t612*t98 2265 t625 = -t82*t479/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t85*t490 - 0.4e1_dp/ & 2266 0.27e2_dp*t88*t494 - 0.32e2_dp/0.243e3_dp*t91*t498 - 0.80e2_dp/0.729e3_dp & 2267 *t94*t502 2268 t627 = t81*t231*t625 2269 t629 = t600*t80*t129 2270 t630 = t614*t118 2271 t632 = t117*R*t625 2272 t633 = t263*t632 2273 t649 = 0.3e1_dp*t272*t273*t600 + 0.3e1_dp*t278*t112*t13*t612 & 2274 - 0.3e1_dp*t285*t112*t13*t625 + t109*(-t601 - t614 + t627) & 2275 *t110*t273 2276 t651 = t269*R*t649 2277 t653 = t99*t651/0.3e1_dp 2278 t655 = (t601 + t614 - t627 - t629 - t630 + t633 + t653)*t121 2279 t658 = t121*t599 2280 t660 = t612*t98 2281 t662 = t312*t625 2282 t674 = (-t601 - t614 + t627 - t629 - t630 + t633 + t653)*t132 2283 t677 = t132*t599 2284 t690 = -0.2e1_dp*t655 + t655*t99 + t658*t308 + t123*t660 - t123 & 2285 *t662 - t655*t81*t129 - t658*t307*t129 - t123*t612* & 2286 t129 + t321*t632 + t125*t651/0.3e1_dp + 0.2e1_dp*t674 + t674*t99 & 2287 + t677*t308 + t134*t660 - t134*t662 + t674*t81*t129 + & 2288 t677*t307*t129 + t134*t612*t129 - t339*t632 - t135*t651 & 2289 /0.3e1_dp 2290 e_laplace_rho = e_laplace_rho + (rho*t690*t117/0.8e1_dp - t139*t269*t649/0.24e2_dp)*sx 2291 END IF 2292 2293 END SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_gt_b 2294 2295! ************************************************************************************************** 2296!> \brief Truncated long range part for y > 0 and R <= b 2297!> \param rho ... 2298!> \param ndrho ... 2299!> \param tau ... 2300!> \param laplace_rho ... 2301!> \param e_0 ... 2302!> \param e_rho ... 2303!> \param e_ndrho ... 2304!> \param e_tau ... 2305!> \param e_laplace_rho ... 2306!> \param sx ... 2307!> \param R ... 2308!> \param gamma ... 2309!> \param grad_deriv ... 2310!> \par History 2311!> 12.2008 created [mguidon] 2312!> \author mguidon 2313! ************************************************************************************************** 2314 SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_lte_b(rho, ndrho, tau, laplace_rho, e_0, & 2315 e_rho, e_ndrho, e_tau, e_laplace_rho, & 2316 sx, R, gamma, grad_deriv) 2317 REAL(dp), INTENT(IN) :: rho, ndrho, tau, laplace_rho 2318 REAL(dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_tau, e_laplace_rho 2319 REAL(dp), INTENT(IN) :: sx, R, gamma 2320 INTEGER, INTENT(IN) :: grad_deriv 2321 2322 CHARACTER(len=*), PARAMETER :: routineN = 'x_br_lsd_y_gt_0_cutoff_R_lte_b', & 2323 routineP = moduleN//':'//routineN 2324 2325 REAL(KIND=dp) :: t1, t100, t101, t102, t103, t104, t105, t106, t108, t109, t110, t111, t112, & 2326 t114, t115, t116, t117, t118, t119, t12, t121, t123, t124, t125, t126, t128, t13, t130, & 2327 t131, t133, t134, t136, t138, t144, t151, t154, t157, t158, t161, t163, t164, t175, t178, & 2328 t179, t180, t182, t184, t185, t187, t19, t190, t192, t193, t2, t20, t201, t205, t206, & 2329 t209, t210, t213, t215, t216, t220, t221, t224, t227, t228, t229, t23, t232, t234, t235, & 2330 t237, t238, t24, t25, t26, t264, t267, t27, t270, t271, t274, t275, t28, t280, t285, & 2331 t287, t29, t292, t294, t3, t30, t302, t305, t308, t31, t310, t315 2332 REAL(KIND=dp) :: t318, t319, t330, t333, t336, t337, t340, t348, t35, t353, t356, t357, t36, & 2333 t365, t366, t368, t37, t371, t375, t379, t383, t387, t390, t392, t393, t4, t41, t410, & 2334 t413, t42, t426, t428, t43, t433, t436, t439, t44, t445, t45, t46, t461, t462, t465, & 2335 t479, t480, t485, t486, t488, t49, t492, t493, t496, t497, t5, t500, t501, t504, t505, & 2336 t508, t510, t511, t52, t523, t526, t53, t539, t54, t541, t546, t549, t55, t552, t558, & 2337 t57, t574, t575, t578, t58, t597, t598, t6, t600, t61, t612, t614, t615, t62, t627, t63, & 2338 t630, t64, t643, t645, t65, t650, t653, t656, t66, t662, t67, t678 2339 REAL(KIND=dp) :: t679, t682, t7, t71, t72, t73, t74, t76, t77, t8, t80, t81, t82, t85, t88, & 2340 t9, t91, t94, t97, t98, t99 2341 2342 IF (grad_deriv >= 0) THEN 2343 t1 = 0.1e1_dp/br_BB 2344 t2 = pi**(0.1e1_dp/0.3e1_dp) 2345 t3 = t2**2 2346 t4 = 0.1e1_dp/t3 2347 t5 = t1*t4 2348 t6 = rho**(0.1e1_dp/0.3e1_dp) 2349 t7 = t6**2 2350 t8 = t7*rho 2351 t9 = 0.1e1_dp/t8 2352 t12 = ndrho**2 2353 t13 = 0.1e1_dp/rho 2354 t19 = laplace_rho/0.6e1_dp - gamma*(REAL(2*tau, KIND=dp) - t12*t13/0.4e1_dp)/0.3e1_dp 2355 t20 = t9*t19 2356 t23 = br_BB**2 2357 t24 = t2*pi 2358 t25 = t23*t24 2359 t26 = rho**2 2360 t27 = t26*rho 2361 t28 = t6*t27 2362 t29 = t19**2 2363 t30 = 0.1e1_dp/t29 2364 t31 = t28*t30 2365 t35 = SQRT(0.10e1_dp + 0.4e1_dp/0.9e1_dp*t25*t31) 2366 t36 = t35*t1 2367 t37 = t4*t9 2368 t41 = 0.1500000000000000e1_dp*t5*t20 + 0.3e1_dp/0.2e1_dp*t36*t37* & 2369 t19 2370 t42 = LOG(t41) 2371 t43 = t42 + 0.2e1_dp 2372 t44 = br_d1*t3 2373 t45 = 0.1e1_dp/t19 2374 t46 = t8*t45 2375 t49 = br_d2*t24 2376 t52 = pi**2 2377 t53 = br_d3*t52 2378 t54 = t26**2 2379 t55 = t54*rho 2380 t57 = 0.1e1_dp/t29/t19 2381 t58 = t55*t57 2382 t61 = t3*t52 2383 t62 = br_d4*t61 2384 t63 = t54*t26 2385 t64 = t7*t63 2386 t65 = t29**2 2387 t66 = 0.1e1_dp/t65 2388 t67 = t64*t66 2389 t71 = t2*t52*pi 2390 t72 = br_d5*t71 2391 t73 = t54**2 2392 t74 = t6*t73 2393 t76 = 0.1e1_dp/t65/t19 2394 t77 = t74*t76 2395 t80 = br_d0 + 0.2e1_dp/0.3e1_dp*t44*t46 + 0.4e1_dp/0.9e1_dp*t49*t31 & 2396 + 0.8e1_dp/0.27e2_dp*t53*t58 + 0.16e2_dp/0.81e2_dp*t62*t67 + 0.32e2_dp & 2397 /0.243e3_dp*t72*t77 2398 t81 = t43*t80 2399 t82 = br_e1*t3 2400 t85 = br_e2*t24 2401 t88 = br_e3*t52 2402 t91 = br_e4*t61 2403 t94 = br_e5*t71 2404 t97 = br_e0 + 0.2e1_dp/0.3e1_dp*t82*t46 + 0.4e1_dp/0.9e1_dp*t85*t31 & 2405 + 0.8e1_dp/0.27e2_dp*t88*t58 + 0.16e2_dp/0.81e2_dp*t91*t67 + 0.32e2_dp & 2406 /0.243e3_dp*t94*t77 2407 t98 = 0.1e1_dp/t97 2408 t99 = t81*t98 2409 t100 = 8**(0.1e1_dp/0.3e1_dp) 2410 t101 = t43**2 2411 t102 = t101*t43 2412 t103 = t80**2 2413 t104 = t103*t80 2414 t105 = t102*t104 2415 t106 = t97**2 2416 t108 = 0.1e1_dp/t106/t97 2417 t109 = t105*t108 2418 t110 = EXP(-t99) 2419 t111 = 0.1e1_dp/0.3141592654e1_dp 2420 t112 = t110*t111 2421 t114 = t109*t112*t13 2422 t115 = t114**(0.1e1_dp/0.3e1_dp) 2423 t116 = 0.1e1_dp/t115 2424 t117 = REAL(t100, KIND=dp)*t116 2425 t118 = t117*R 2426 t119 = t99*t118 2427 t121 = EXP(0.2e1_dp*t119) 2428 t123 = t121*R 2429 t124 = t123*t43 2430 t125 = t80*t98 2431 t126 = t125*t117 2432 t128 = t121*t43 2433 t130 = t99 + t119 2434 t131 = EXP(t130) 2435 t133 = 0.2e1_dp*t121 - t124*t126 + t128*t125 + 0.2e1_dp + t99 + t119 & 2436 - 0.4e1_dp*t131 2437 t134 = rho*t133 2438 t136 = EXP(-t130) 2439 t138 = t136*REAL(t100, KIND=dp)*t116 2440 e_0 = e_0 + (t134*t138/0.8e1_dp)*sx 2441 END IF 2442 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 2443 t144 = 0.1e1_dp/t7/t26 2444 t151 = 0.1e1_dp/t7/t27*gamma*t12 2445 t154 = 0.1e1_dp/t35 2446 t157 = t6*t26 2447 t158 = t157*t30 2448 t161 = t6*rho 2449 t163 = t57*gamma 2450 t164 = t163*t12 2451 t175 = t36*t4 2452 t178 = -0.2500000000e1_dp*t5*t144*t19 - 0.1250000000e0_dp*t5*t151 & 2453 + 0.3e1_dp/0.4e1_dp*t154*t1*t4*t20*(0.40e2_dp/0.27e2_dp*t25 & 2454 *t158 + 0.2e1_dp/0.27e2_dp*t25*t161*t164) - 0.5e1_dp/0.2e1_dp*t36 & 2455 *t4*t144*t19 - t175*t151/0.8e1_dp 2456 t179 = 0.1e1_dp/t41 2457 t180 = t178*t179 2458 t182 = t98*REAL(t100, KIND=dp) 2459 t184 = t182*t116*R 2460 t185 = t180*t80*t184 2461 t187 = t7*t45 2462 t190 = 0.1e1_dp/t6 2463 t192 = t30*gamma 2464 t193 = t192*t12 2465 t201 = t54*t57 2466 t205 = t66*gamma 2467 t206 = t205*t12 2468 t209 = t7*t55 2469 t210 = t209*t66 2470 t213 = t7*t54 2471 t215 = t76*gamma 2472 t216 = t215*t12 2473 t220 = t6*t54*t27 2474 t221 = t220*t76 2475 t224 = t6*t63 2476 t227 = 0.1e1_dp/t65/t29 2477 t228 = t227*gamma 2478 t229 = t228*t12 2479 t232 = 0.10e2_dp/0.9e1_dp*t44*t187 + t44*t190*t193/0.18e2_dp + 0.40e2_dp & 2480 /0.27e2_dp*t49*t158 + 0.2e1_dp/0.27e2_dp*t49*t161*t164 + & 2481 0.40e2_dp/0.27e2_dp*t53*t201 + 0.2e1_dp/0.27e2_dp*t53*t27*t206 & 2482 + 0.320e3_dp/0.243e3_dp*t62*t210 + 0.16e2_dp/0.243e3_dp*t62*t213* & 2483 t216 + 0.800e3_dp/0.729e3_dp*t72*t221 + 0.40e2_dp/0.729e3_dp*t72* & 2484 t224*t229 2485 t234 = t43*t232*t98 2486 t235 = t234*t118 2487 t237 = 0.1e1_dp/t106 2488 t238 = t81*t237 2489 t264 = 0.10e2_dp/0.9e1_dp*t82*t187 + t82*t190*t193/0.18e2_dp + 0.40e2_dp & 2490 /0.27e2_dp*t85*t158 + 0.2e1_dp/0.27e2_dp*t85*t161*t164 + & 2491 0.40e2_dp/0.27e2_dp*t88*t201 + 0.2e1_dp/0.27e2_dp*t88*t27*t206 & 2492 + 0.320e3_dp/0.243e3_dp*t91*t210 + 0.16e2_dp/0.243e3_dp*t91*t213* & 2493 t216 + 0.800e3_dp/0.729e3_dp*t94*t221 + 0.40e2_dp/0.729e3_dp*t94* & 2494 t224*t229 2495 t267 = t238*t117*R*t264 2496 t270 = 0.1e1_dp/t115/t114 2497 t271 = REAL(t100, KIND=dp)*t270 2498 t274 = t101*t104*t108*t110 2499 t275 = t111*t13 2500 t280 = t102*t103*t108 2501 t285 = t106**2 2502 t287 = t105/t285 2503 t292 = t180*t125 2504 t294 = t81*t237*t264 2505 t302 = 0.3e1_dp*t274*t275*t180 + 0.3e1_dp*t280*t112*t13*t232 & 2506 - 0.3e1_dp*t287*t112*t13*t264 + t109*(-t292 - t234 + t294) & 2507 *t110*t275 - t109*t112/t26 2508 t305 = t99*t271*R*t302 2509 t308 = (0.2e1_dp*t185 + 0.2e1_dp*t235 - 0.2e1_dp*t267 - 0.2e1_dp/0.3e1_dp & 2510 *t305)*t121 2511 t310 = R*t43 2512 t315 = t232*t98 2513 t318 = t123*t81 2514 t319 = t237*REAL(t100, KIND=dp) 2515 t330 = t179*t80*t98 2516 t333 = t80*t237 2517 t336 = t305/0.3e1_dp 2518 t337 = t292 + t234 - t294 + t185 + t235 - t267 - t336 2519 t340 = 0.2e1_dp*t308 - t308*t310*t126 - t123*t180*t126 - t124 & 2520 *t315*t117 + t318*t319*t116*t264 + t318*t182*t270* & 2521 t302/0.3e1_dp + t308*t99 + t121*t178*t330 + t128*t315 - t128 & 2522 *t333*t264 + t292 + t234 - t294 + t185 + t235 - t267 - t336 & 2523 - 0.4e1_dp*t337*t131 2524 t348 = t134*t136 2525 e_rho = e_rho + (t133*t136*t117/0.8e1_dp + rho*t340*t138/0.8e1_dp - t134 & 2526 *t337*t138/0.8e1_dp - t348*t271*t302/0.24e2_dp)*sx 2527 t353 = t144*gamma*ndrho 2528 t356 = t154*br_BB 2529 t357 = t356*t3 2530 t365 = 0.2500000000000000e0_dp*t5*t353 - t357*t7*t30*gamma*ndrho & 2531 /0.9e1_dp + t175*t353/0.4e1_dp 2532 t366 = t365*t179 2533 t368 = t366*t80*t184 2534 t371 = t192*ndrho 2535 t375 = t163*ndrho 2536 t379 = t205*ndrho 2537 t383 = t215*ndrho 2538 t387 = t228*ndrho 2539 t390 = -t44*t7*t371/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t49*t157*t375 & 2540 - 0.4e1_dp/0.27e2_dp*t53*t54*t379 - 0.32e2_dp/0.243e3_dp*t62*t209 & 2541 *t383 - 0.80e2_dp/0.729e3_dp*t72*t220*t387 2542 t392 = t43*t390*t98 2543 t393 = t392*t118 2544 t410 = -t82*t7*t371/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t85*t157*t375 & 2545 - 0.4e1_dp/0.27e2_dp*t88*t54*t379 - 0.32e2_dp/0.243e3_dp*t91*t209 & 2546 *t383 - 0.80e2_dp/0.729e3_dp*t94*t220*t387 2547 t413 = t238*t117*R*t410 2548 t426 = t366*t125 2549 t428 = t81*t237*t410 2550 t433 = 0.3e1_dp*t274*t275*t366 + 0.3e1_dp*t280*t112*t13*t390 & 2551 - 0.3e1_dp*t287*t112*t13*t410 + t109*(-t426 - t392 + t428) & 2552 *t110*t275 2553 t436 = t99*t271*R*t433 2554 t439 = (0.2e1_dp*t368 + 0.2e1_dp*t393 - 0.2e1_dp*t413 - 0.2e1_dp/0.3e1_dp & 2555 *t436)*t121 2556 t445 = t390*t98 2557 t461 = t436/0.3e1_dp 2558 t462 = t426 + t392 - t428 + t368 + t393 - t413 - t461 2559 t465 = 0.2e1_dp*t439 - t439*t310*t126 - t123*t366*t126 - t124 & 2560 *t445*t117 + t318*t319*t116*t410 + t318*t182*t270* & 2561 t433/0.3e1_dp + t439*t99 + t121*t365*t330 + t128*t445 - t128 & 2562 *t333*t410 + t426 + t392 - t428 + t368 + t393 - t413 - t461 & 2563 - 0.4e1_dp*t462*t131 2564 e_ndrho = e_ndrho + (rho*t465*t138/0.8e1_dp - t134*t462*t138/0.8e1_dp - t348 & 2565 *t271*t433/0.24e2_dp)*sx 2566 t479 = t8*t30 2567 t480 = t479*gamma 2568 t485 = -0.1000000000e1_dp*t5*t9*gamma + 0.4e1_dp/0.9e1_dp*t357*t480 & 2569 - t36*t37*gamma 2570 t486 = t485*t179 2571 t488 = t486*t80*t184 2572 t492 = t28*t57 2573 t493 = t492*gamma 2574 t496 = t55*t66 2575 t497 = t496*gamma 2576 t500 = t64*t76 2577 t501 = t500*gamma 2578 t504 = t74*t227 2579 t505 = t504*gamma 2580 t508 = 0.4e1_dp/0.9e1_dp*t44*t480 + 0.16e2_dp/0.27e2_dp*t49*t493 + & 2581 0.16e2_dp/0.27e2_dp*t53*t497 + 0.128e3_dp/0.243e3_dp*t62*t501 + 0.320e3_dp & 2582 /0.729e3_dp*t72*t505 2583 t510 = t43*t508*t98 2584 t511 = t510*t118 2585 t523 = 0.4e1_dp/0.9e1_dp*t82*t480 + 0.16e2_dp/0.27e2_dp*t85*t493 + & 2586 0.16e2_dp/0.27e2_dp*t88*t497 + 0.128e3_dp/0.243e3_dp*t91*t501 + 0.320e3_dp & 2587 /0.729e3_dp*t94*t505 2588 t526 = t238*t117*R*t523 2589 t539 = t486*t125 2590 t541 = t81*t237*t523 2591 t546 = 0.3e1_dp*t274*t275*t486 + 0.3e1_dp*t280*t112*t13*t508 & 2592 - 0.3e1_dp*t287*t112*t13*t523 + t109*(-t539 - t510 + t541) & 2593 *t110*t275 2594 t549 = t99*t271*R*t546 2595 t552 = (0.2e1_dp*t488 + 0.2e1_dp*t511 - 0.2e1_dp*t526 - 0.2e1_dp/0.3e1_dp & 2596 *t549)*t121 2597 t558 = t508*t98 2598 t574 = t549/0.3e1_dp 2599 t575 = t539 + t510 - t541 + t488 + t511 - t526 - t574 2600 t578 = 0.2e1_dp*t552 - t552*t310*t126 - t123*t486*t126 - t124 & 2601 *t558*t117 + t318*t319*t116*t523 + t318*t182*t270* & 2602 t546/0.3e1_dp + t552*t99 + t121*t485*t330 + t128*t558 - t128 & 2603 *t333*t523 + t539 + t510 - t541 + t488 + t511 - t526 - t574 & 2604 - 0.4e1_dp*t575*t131 2605 e_tau = e_tau + (rho*t578*t138/0.8e1_dp - t134*t575*t138/0.8e1_dp - t348 & 2606 *t271*t546/0.24e2_dp)*sx 2607 t597 = 0.2500000000000000e0_dp*t5*t9 - t356*t3*t8*t30/0.9e1_dp & 2608 + t36*t37/0.4e1_dp 2609 t598 = t597*t179 2610 t600 = t598*t80*t184 2611 t612 = -t44*t479/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t49*t492 - 0.4e1_dp/ & 2612 0.27e2_dp*t53*t496 - 0.32e2_dp/0.243e3_dp*t62*t500 - 0.80e2_dp/0.729e3_dp & 2613 *t72*t504 2614 t614 = t43*t612*t98 2615 t615 = t614*t118 2616 t627 = -t82*t479/0.9e1_dp - 0.4e1_dp/0.27e2_dp*t85*t492 - 0.4e1_dp/ & 2617 0.27e2_dp*t88*t496 - 0.32e2_dp/0.243e3_dp*t91*t500 - 0.80e2_dp/0.729e3_dp & 2618 *t94*t504 2619 t630 = t238*t117*R*t627 2620 t643 = t598*t125 2621 t645 = t81*t237*t627 2622 t650 = 0.3e1_dp*t274*t275*t598 + 0.3e1_dp*t280*t112*t13*t612 & 2623 - 0.3e1_dp*t287*t112*t13*t627 + t109*(-t643 - t614 + t645) & 2624 *t110*t275 2625 t653 = t99*t271*R*t650 2626 t656 = (0.2e1_dp*t600 + 0.2e1_dp*t615 - 0.2e1_dp*t630 - 0.2e1_dp/0.3e1_dp & 2627 *t653)*t121 2628 t662 = t612*t98 2629 t678 = t653/0.3e1_dp 2630 t679 = t643 + t614 - t645 + t600 + t615 - t630 - t678 2631 t682 = 0.2e1_dp*t656 - t656*t310*t126 - t123*t598*t126 - t124 & 2632 *t662*t117 + t318*t319*t116*t627 + t318*t182*t270* & 2633 t650/0.3e1_dp + t656*t99 + t121*t597*t330 + t128*t662 - t128 & 2634 *t333*t627 + t643 + t614 - t645 + t600 + t615 - t630 - t678 & 2635 - 0.4e1_dp*t679*t131 2636 e_laplace_rho = e_laplace_rho + (rho*t682*t138/0.8e1_dp - t134*t679*t138/0.8e1_dp - t348 & 2637 *t271*t650/0.24e2_dp)*sx 2638 END IF 2639 2640 END SUBROUTINE x_br_lsd_y_gt_0_cutoff_R_lte_b 2641 2642END MODULE xc_xbecke_roussel 2643