1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Set of procedures to modify the density matrix, in order to collocate 8!> on the real space grid specific functions of the primitives 9!> \author MI 06.2006 10! ************************************************************************************************** 11MODULE grid_modify_pab_block 12 13 USE kinds, ONLY: dp 14 USE orbital_pointers, ONLY: coset 15#include "../base/base_uses.f90" 16 17 IMPLICIT NONE 18 19 PRIVATE 20 21 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'grid_modify_pab_block' 22 23 PUBLIC :: prepare_dadb, prepare_adb_m_dab, prepare_dab_p_adb, prepare_ardb_m_darb 24 PUBLIC :: prepare_diadib, prepare_dijadijb, prepare_diiadiib 25 26CONTAINS 27 28! ************************************************************************************************** 29!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 30!> is equivalent to mapping pab with 0.5 * (nabla pgf_a) . (nabla pgf_b) 31!> (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) 32!> \param pab_local ... 33!> \param pab ... 34!> \param lxa ... 35!> \param lya ... 36!> \param lza ... 37!> \param lxb ... 38!> \param lyb ... 39!> \param lzb ... 40!> \param o1 ... 41!> \param o2 ... 42!> \param zeta ... 43!> \param zetb ... 44! ************************************************************************************************** 45 SUBROUTINE prepare_dadb(pab_local, pab, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 46 47 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 48 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 49 INTEGER, INTENT(IN) :: lxa, lya, lza, lxb, lyb, lzb, o1, o2 50 REAL(dp), INTENT(IN) :: zeta, zetb 51 52 INTEGER :: ico, ico_l, jco, jco_l 53 54! this element of pab results in 12 elements of pab_local 55 56 ico = coset(lxa, lya, lza) 57 jco = coset(lxb, lyb, lzb) 58 ! x (all safe if lxa = 0, as the spurious added terms have zero prefactor) 59 60 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 61 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 62 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxa*lxb*pab(o1 + ico, o2 + jco) 63 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 64 jco_l = coset((lxb + 1), lyb, lzb) 65 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lxa*zetb*pab(o1 + ico, o2 + jco) 66 ico_l = coset((lxa + 1), lya, lza) 67 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 68 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lxb*pab(o1 + ico, o2 + jco) 69 ico_l = coset((lxa + 1), lya, lza) 70 jco_l = coset((lxb + 1), lyb, lzb) 71 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 72 73 ! y 74 75 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 76 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 77 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lya*lyb*pab(o1 + ico, o2 + jco) 78 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 79 jco_l = coset(lxb, (lyb + 1), lzb) 80 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lya*zetb*pab(o1 + ico, o2 + jco) 81 ico_l = coset(lxa, (lya + 1), lza) 82 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 83 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lyb*pab(o1 + ico, o2 + jco) 84 ico_l = coset(lxa, (lya + 1), lza) 85 jco_l = coset(lxb, (lyb + 1), lzb) 86 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 87 88 ! z 89 90 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 91 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 92 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lza*lzb*pab(o1 + ico, o2 + jco) 93 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 94 jco_l = coset(lxb, lyb, (lzb + 1)) 95 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lza*zetb*pab(o1 + ico, o2 + jco) 96 ico_l = coset(lxa, lya, (lza + 1)) 97 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 98 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lzb*pab(o1 + ico, o2 + jco) 99 ico_l = coset(lxa, lya, (lza + 1)) 100 jco_l = coset(lxb, lyb, (lzb + 1)) 101 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 102 103 END SUBROUTINE prepare_dadb 104 105! ************************************************************************************************** 106!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 107!> is equivalent to mapping pab with (ddi pgf_a) . (ddi pgf_b) 108!> (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) 109!> \param pab_local ... 110!> \param pab ... 111!> \param ider ... 112!> \param lxa ... 113!> \param lya ... 114!> \param lza ... 115!> \param lxb ... 116!> \param lyb ... 117!> \param lzb ... 118!> \param o1 ... 119!> \param o2 ... 120!> \param zeta ... 121!> \param zetb ... 122! ************************************************************************************************** 123 SUBROUTINE prepare_diadib(pab_local, pab, ider, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 124 125 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 126 REAL(dp), DIMENSION(:, :), INTENT(in) :: pab 127 INTEGER, INTENT(in) :: ider, lxa, lya, lza, lxb, lyb, lzb, o1, & 128 o2 129 REAL(dp), INTENT(in) :: zeta, zetb 130 131 INTEGER :: ico, ico_l, jco, jco_l 132 133! this element of pab results in 12 elements of pab_local 134 135 ico = coset(lxa, lya, lza) 136 jco = coset(lxb, lyb, lzb) 137 IF (ider == 1) THEN 138 ! x (all safe if lxa = 0, as the spurious added terms have zero prefactor) 139 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 140 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 141 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxa*lxb*pab(o1 + ico, o2 + jco) 142 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 143 jco_l = coset((lxb + 1), lyb, lzb) 144 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lxa*zetb*pab(o1 + ico, o2 + jco) 145 ico_l = coset((lxa + 1), lya, lza) 146 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 147 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lxb*pab(o1 + ico, o2 + jco) 148 ico_l = coset((lxa + 1), lya, lza) 149 jco_l = coset((lxb + 1), lyb, lzb) 150 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 151 152 ELSEIF (ider == 2) THEN 153 ! y 154 155 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 156 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 157 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lya*lyb*pab(o1 + ico, o2 + jco) 158 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 159 jco_l = coset(lxb, (lyb + 1), lzb) 160 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lya*zetb*pab(o1 + ico, o2 + jco) 161 ico_l = coset(lxa, (lya + 1), lza) 162 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 163 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lyb*pab(o1 + ico, o2 + jco) 164 ico_l = coset(lxa, (lya + 1), lza) 165 jco_l = coset(lxb, (lyb + 1), lzb) 166 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 167 ELSEIF (ider == 3) THEN 168 ! z 169 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 170 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 171 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lza*lzb*pab(o1 + ico, o2 + jco) 172 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 173 jco_l = coset(lxb, lyb, (lzb + 1)) 174 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lza*zetb*pab(o1 + ico, o2 + jco) 175 ico_l = coset(lxa, lya, (lza + 1)) 176 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 177 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lzb*pab(o1 + ico, o2 + jco) 178 ico_l = coset(lxa, lya, (lza + 1)) 179 jco_l = coset(lxb, lyb, (lzb + 1)) 180 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 181 END IF 182 183 END SUBROUTINE prepare_diadib 184 185! ************************************************************************************************** 186!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 187!> is equivalent to mapping pab with (ddidj pgf_a) . (ddidj pgf_b) 188!> (ddxdy pgf_a ) (ddxdy pgf_b) = 189!> \param pab_local ... 190!> \param pab ... 191!> \param ider1 ... 192!> \param ider2 ... 193!> \param lxa ... 194!> \param lya ... 195!> \param lza ... 196!> \param lxb ... 197!> \param lyb ... 198!> \param lzb ... 199!> \param o1 ... 200!> \param o2 ... 201!> \param zeta ... 202!> \param zetb ... 203! ************************************************************************************************** 204 SUBROUTINE prepare_dijadijb(pab_local, pab, ider1, ider2, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 205 206 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 207 REAL(dp), DIMENSION(:, :), INTENT(in) :: pab 208 INTEGER, INTENT(in) :: ider1, ider2, lxa, lya, lza, lxb, lyb, & 209 lzb, o1, o2 210 REAL(dp), INTENT(in) :: zeta, zetb 211 212 INTEGER :: ico, ico_l, jco 213 REAL(dp) :: func_a 214 215! this element of pab results in 12 elements of pab_local 216 217 ico = coset(lxa, lya, lza) 218 jco = coset(lxb, lyb, lzb) 219 220 IF ((ider1 == 1 .AND. ider2 == 2) .OR. (ider1 == 2 .AND. ider2 == 1)) THEN ! xy 221 ico_l = coset(MAX(lxa - 1, 0), MAX(lya - 1, 0), lza) 222 func_a = lxa*lya*pab(o1 + ico, o2 + jco) 223 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 224 ico_l = coset(lxa + 1, MAX(lya - 1, 0), lza) 225 func_a = -2.0_dp*zeta*lya*pab(o1 + ico, o2 + jco) 226 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 227 ico_l = coset(MAX(lxa - 1, 0), lya + 1, lza) 228 func_a = -2.0_dp*zeta*lxa*pab(o1 + ico, o2 + jco) 229 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 230 ico_l = coset(lxa + 1, lya + 1, lza) 231 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 232 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 233 ELSEIF ((ider1 == 2 .AND. ider2 == 3) .OR. (ider1 == 3 .AND. ider2 == 2)) THEN ! yz 234 ico_l = coset(lxa, MAX(lya - 1, 0), MAX(lza - 1, 0)) 235 func_a = lya*lza*pab(o1 + ico, o2 + jco) 236 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 237 ico_l = coset(lxa, lya + 1, MAX(lza - 1, 0)) 238 func_a = -2.0_dp*zeta*lza*pab(o1 + ico, o2 + jco) 239 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 240 ico_l = coset(lxa, MAX(lya - 1, 0), lza + 1) 241 func_a = -2.0_dp*zeta*lya*pab(o1 + ico, o2 + jco) 242 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 243 ico_l = coset(lxa, lya + 1, lza + 1) 244 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 245 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 246 ELSEIF ((ider1 == 3 .AND. ider2 == 1) .OR. (ider1 == 1 .AND. ider2 == 3)) THEN ! zx 247 ico_l = coset(MAX(lxa - 1, 0), lya, MAX(lza - 1, 0)) 248 func_a = lza*lxa*pab(o1 + ico, o2 + jco) 249 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 250 ico_l = coset(MAX(lxa - 1, 0), lya, lza + 1) 251 func_a = -2.0_dp*zeta*lxa*pab(o1 + ico, o2 + jco) 252 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 253 ico_l = coset(lxa + 1, lya, MAX(lza - 1, 0)) 254 func_a = -2.0_dp*zeta*lza*pab(o1 + ico, o2 + jco) 255 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 256 ico_l = coset(lxa + 1, lya, lza + 1) 257 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 258 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 259 END IF 260 261 END SUBROUTINE prepare_dijadijb 262 263! ************************************************************************************************** 264!> \brief ... 265!> \param pab_local ... 266!> \param func_a ... 267!> \param ico_l ... 268!> \param lx ... 269!> \param ly ... 270!> \param lz ... 271!> \param zet ... 272!> \param idir ... 273! ************************************************************************************************** 274 SUBROUTINE oneterm_dijdij(pab_local, func_a, ico_l, lx, ly, lz, zet, idir) 275 276 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 277 REAL(dp), INTENT(in) :: func_a 278 INTEGER, INTENT(in) :: ico_l, lx, ly, lz 279 REAL(dp), INTENT(in) :: zet 280 INTEGER, INTENT(in) :: idir 281 282 INTEGER :: jco_l, l1, l2 283 284 IF (idir == 1) THEN 285 l1 = lx 286 l2 = ly 287 jco_l = coset(MAX(lx - 1, 0), MAX(ly - 1, 0), lz) 288 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*l2*func_a 289 jco_l = coset(lx + 1, MAX(ly - 1, 0), lz) 290 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l2*func_a 291 jco_l = coset(MAX(lx - 1, 0), ly + 1, lz) 292 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l1*func_a 293 jco_l = coset(lx + 1, ly + 1, lz) 294 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 295 ELSEIF (idir == 2) THEN 296 l1 = ly 297 l2 = lz 298 jco_l = coset(lx, MAX(ly - 1, 0), MAX(lz - 1, 0)) 299 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*l2*func_a 300 jco_l = coset(lx, ly + 1, MAX(lz - 1, 0)) 301 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l2*func_a 302 jco_l = coset(lx, MAX(ly - 1, 0), lz + 1) 303 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l1*func_a 304 jco_l = coset(lx, ly + 1, lz + 1) 305 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 306 ELSEIF (idir == 3) THEN 307 l1 = lz 308 l2 = lx 309 jco_l = coset(MAX(lx - 1, 0), ly, MAX(lz - 1, 0)) 310 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*l2*func_a 311 jco_l = coset(MAX(lx - 1, 0), ly, lz + 1) 312 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l2*func_a 313 jco_l = coset(lx + 1, ly, MAX(lz - 1, 0)) 314 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l1*func_a 315 jco_l = coset(lx + 1, ly, lz + 1) 316 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 317 318 END IF 319 END SUBROUTINE oneterm_dijdij 320 321! 322! ************************************************************************************************** 323!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 324!> is equivalent to mapping pab with (ddidj pgf_a) . (ddidj pgf_b) 325!> (ddxdx pgf_a ) (ddxdx pgf_b) = 326!> \param pab_local ... 327!> \param pab ... 328!> \param ider ... 329!> \param lxa ... 330!> \param lya ... 331!> \param lza ... 332!> \param lxb ... 333!> \param lyb ... 334!> \param lzb ... 335!> \param o1 ... 336!> \param o2 ... 337!> \param zeta ... 338!> \param zetb ... 339! ************************************************************************************************** 340 SUBROUTINE prepare_diiadiib(pab_local, pab, ider, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 341 342 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 343 REAL(dp), DIMENSION(:, :), INTENT(in) :: pab 344 INTEGER, INTENT(in) :: ider, lxa, lya, lza, lxb, lyb, lzb, o1, & 345 o2 346 REAL(dp), INTENT(in) :: zeta, zetb 347 348 INTEGER :: ico, ico_l, jco 349 REAL(dp) :: func_a 350 351! this element of pab results in 9 elements of pab_local 352 353 ico = coset(lxa, lya, lza) 354 jco = coset(lxb, lyb, lzb) 355 356 IF (ider == 1) THEN ! x 357 ico_l = coset(MAX(lxa - 2, 0), lya, lza) 358 func_a = lxa*(lxa - 1)*pab(o1 + ico, o2 + jco) 359 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 360 ico_l = coset(lxa, lya, lza) 361 func_a = -2.0_dp*zeta*(2*lxa + 1)*pab(o1 + ico, o2 + jco) 362 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 363 ico_l = coset(lxa + 2, lya, lza) 364 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 365 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 366 ELSEIF (ider == 2) THEN ! y 367 ico_l = coset(lxa, MAX(lya - 2, 0), lza) 368 func_a = lya*(lya - 1)*pab(o1 + ico, o2 + jco) 369 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 370 ico_l = coset(lxa, lya, lza) 371 func_a = -2.0_dp*zeta*(2*lya + 1)*pab(o1 + ico, o2 + jco) 372 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 373 ico_l = coset(lxa, lya + 2, lza) 374 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 375 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 376 ELSEIF (ider == 3) THEN ! z 377 ico_l = coset(lxa, lya, MAX(lza - 2, 0)) 378 func_a = lza*(lza - 1)*pab(o1 + ico, o2 + jco) 379 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 380 ico_l = coset(lxa, lya, lza) 381 func_a = -2.0_dp*zeta*(2*lza + 1)*pab(o1 + ico, o2 + jco) 382 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 383 ico_l = coset(lxa, lya, lza + 2) 384 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 385 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 386 END IF 387 388 END SUBROUTINE prepare_diiadiib 389 390! ************************************************************************************************** 391!> \brief ... 392!> \param pab_local ... 393!> \param func_a ... 394!> \param ico_l ... 395!> \param lx ... 396!> \param ly ... 397!> \param lz ... 398!> \param zet ... 399!> \param idir ... 400! ************************************************************************************************** 401 SUBROUTINE oneterm_diidii(pab_local, func_a, ico_l, lx, ly, lz, zet, idir) 402 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 403 REAL(dp), INTENT(in) :: func_a 404 INTEGER, INTENT(in) :: ico_l, lx, ly, lz 405 REAL(dp), INTENT(in) :: zet 406 INTEGER, INTENT(in) :: idir 407 408 INTEGER :: jco_l, l1 409 410 IF (idir == 1) THEN 411 l1 = lx 412 jco_l = coset(MAX(lx - 2, 0), ly, lz) 413 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*(l1 - 1)*func_a 414 jco_l = coset(lx, ly, lz) 415 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*(2*l1 + 1)*func_a 416 jco_l = coset(lx + 2, ly, lz) 417 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 418 ELSEIF (idir == 2) THEN 419 l1 = ly 420 jco_l = coset(lx, MAX(ly - 2, 0), lz) 421 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*(l1 - 1)*func_a 422 jco_l = coset(lx, ly, lz) 423 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*(2*l1 + 1)*func_a 424 jco_l = coset(lx, ly + 2, lz) 425 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 426 ELSEIF (idir == 3) THEN 427 l1 = lz 428 jco_l = coset(lx, ly, MAX(lz - 2, 0)) 429 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*(l1 - 1)*func_a 430 jco_l = coset(lx, ly, lz) 431 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*(2*l1 + 1)*func_a 432 jco_l = coset(lx, ly, lz + 2) 433 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 434 END IF 435 436 END SUBROUTINE oneterm_diidii 437 438! ************************************************************************************************** 439!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 440!> is equivalent to mapping pab with pgf_a (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) pgf_b 441!> ( pgf_a ) (ddx pgf_b) - (ddx pgf_a)( pgf_b ) = 442!> pgf_a *(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) - (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b 443!> \param pab_local ... 444!> \param pab ... 445!> \param idir ... 446!> \param lxa ... 447!> \param lya ... 448!> \param lza ... 449!> \param lxb ... 450!> \param lyb ... 451!> \param lzb ... 452!> \param o1 ... 453!> \param o2 ... 454!> \param zeta ... 455!> \param zetb ... 456! ************************************************************************************************** 457 SUBROUTINE prepare_adb_m_dab(pab_local, pab, idir, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 458 459 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 460 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 461 INTEGER, INTENT(IN) :: idir, lxa, lya, lza, lxb, lyb, lzb, o1, & 462 o2 463 REAL(dp), INTENT(IN) :: zeta, zetb 464 465 INTEGER :: ico, ico_l, jco, jco_l 466 467! this element of pab results in 4 elements of pab_local 468 469 ico = coset(lxa, lya, lza) 470 jco = coset(lxb, lyb, lzb) 471 472 IF (idir == 1) THEN ! x 473 ico_l = coset(lxa, lya, lza) 474 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 475 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 476 ico_l = coset(lxa, lya, lza) 477 jco_l = coset((lxb + 1), lyb, lzb) 478 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 479 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 480 jco_l = coset(lxb, lyb, lzb) 481 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lxa*pab(o1 + ico, o2 + jco) 482 ico_l = coset((lxa + 1), lya, lza) 483 jco_l = coset(lxb, lyb, lzb) 484 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 485 ELSEIF (idir == 2) THEN ! y 486 ico_l = coset(lxa, lya, lza) 487 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 488 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 489 ico_l = coset(lxa, lya, lza) 490 jco_l = coset(lxb, (lyb + 1), lzb) 491 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 492 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 493 jco_l = coset(lxb, lyb, lzb) 494 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lya*pab(o1 + ico, o2 + jco) 495 ico_l = coset(lxa, (lya + 1), lza) 496 jco_l = coset(lxb, lyb, lzb) 497 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 498 ELSE ! z 499 ico_l = coset(lxa, lya, lza) 500 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 501 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 502 ico_l = coset(lxa, lya, lza) 503 jco_l = coset(lxb, lyb, (lzb + 1)) 504 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 505 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 506 jco_l = coset(lxb, lyb, lzb) 507 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lza*pab(o1 + ico, o2 + jco) 508 ico_l = coset(lxa, lya, (lza + 1)) 509 jco_l = coset(lxb, lyb, lzb) 510 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 511 END IF 512 513 END SUBROUTINE prepare_adb_m_dab 514 515! ************************************************************************************************** 516!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 517!> is equivalent to mapping pab with pgf_a (nabla_{idir} pgf_b) + (nabla_{idir} pgf_a) pgf_b 518!> ( pgf_a ) (ddx pgf_b) + (ddx pgf_a)( pgf_b ) = 519!> pgf_a *(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) + (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b 520!> \param pab_local ... 521!> \param pab ... 522!> \param idir ... 523!> \param lxa ... 524!> \param lya ... 525!> \param lza ... 526!> \param lxb ... 527!> \param lyb ... 528!> \param lzb ... 529!> \param o1 ... 530!> \param o2 ... 531!> \param zeta ... 532!> \param zetb ... 533! ************************************************************************************************** 534 SUBROUTINE prepare_dab_p_adb(pab_local, pab, idir, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 535 536 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 537 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 538 INTEGER, INTENT(IN) :: idir, lxa, lya, lza, lxb, lyb, lzb, o1, & 539 o2 540 REAL(dp), INTENT(IN) :: zeta, zetb 541 542 INTEGER :: ico, ico_l, jco, jco_l 543 544! this element of pab results in 4 elements of pab_local 545 546 ico = coset(lxa, lya, lza) 547 jco = coset(lxb, lyb, lzb) 548 549 IF (idir == 1) THEN ! x 550 ico_l = coset(lxa, lya, lza) 551 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 552 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 553 ico_l = coset(lxa, lya, lza) 554 jco_l = coset((lxb + 1), lyb, lzb) 555 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 556 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 557 jco_l = coset(lxb, lyb, lzb) 558 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxa*pab(o1 + ico, o2 + jco) 559 ico_l = coset((lxa + 1), lya, lza) 560 jco_l = coset(lxb, lyb, lzb) 561 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 562 ELSEIF (idir == 2) THEN ! y 563 ico_l = coset(lxa, lya, lza) 564 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 565 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 566 ico_l = coset(lxa, lya, lza) 567 jco_l = coset(lxb, (lyb + 1), lzb) 568 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 569 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 570 jco_l = coset(lxb, lyb, lzb) 571 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lya*pab(o1 + ico, o2 + jco) 572 ico_l = coset(lxa, (lya + 1), lza) 573 jco_l = coset(lxb, lyb, lzb) 574 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 575 ELSE ! z 576 ico_l = coset(lxa, lya, lza) 577 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 578 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 579 ico_l = coset(lxa, lya, lza) 580 jco_l = coset(lxb, lyb, (lzb + 1)) 581 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 582 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 583 jco_l = coset(lxb, lyb, lzb) 584 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lza*pab(o1 + ico, o2 + jco) 585 ico_l = coset(lxa, lya, (lza + 1)) 586 jco_l = coset(lxb, lyb, lzb) 587 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 588 END IF 589 590 END SUBROUTINE prepare_dab_p_adb 591 592! ************************************************************************************************** 593!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 594!> pgf_a (r-Rb)_{ir} (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) (r-Rb)_{ir} pgf_b 595!> ( pgf_a )(r-Rb)_{ir} (ddx pgf_b) - (ddx pgf_a) (r-Rb)_{ir} ( pgf_b ) = 596!> pgf_a *(lbx pgf_{b-1x+1ir} - 2*zetb*pgf_{b+1x+1ir}) - 597!> (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_{b+1ir} 598!> \param pab_local ... 599!> \param pab ... 600!> \param idir ... 601!> \param ir ... 602!> \param lxa ... 603!> \param lya ... 604!> \param lza ... 605!> \param lxb ... 606!> \param lyb ... 607!> \param lzb ... 608!> \param o1 ... 609!> \param o2 ... 610!> \param zeta ... 611!> \param zetb ... 612! ************************************************************************************************** 613 SUBROUTINE prepare_ardb_m_darb(pab_local, pab, idir, ir, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 614 615 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 616 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 617 INTEGER, INTENT(IN) :: idir, ir, lxa, lya, lza, lxb, lyb, lzb, & 618 o1, o2 619 REAL(dp), INTENT(IN) :: zeta, zetb 620 621 INTEGER :: ico, ico_l, jco, jco_l 622 623! this element of pab results in 4 elements of pab_local 624 625 ico = coset(lxa, lya, lza) 626 jco = coset(lxb, lyb, lzb) 627 628 IF (idir == 1 .AND. ir == 1) THEN 629 630 ico_l = coset(lxa, lya, lza) 631 jco_l = coset(lxb, lyb, lzb) 632 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 633 634 ico_l = coset(lxa, lya, lza) 635 jco_l = coset((lxb + 2), lyb, lzb) 636 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 637 638 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 639 jco_l = coset((lxb + 1), lyb, lzb) 640 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lxa*pab(o1 + ico, o2 + jco) 641 642 ico_l = coset((lxa + 1), lya, lza) 643 jco_l = coset((lxb + 1), lyb, lzb) 644 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 645 646 ELSEIF (idir == 1 .AND. ir == 2) THEN 647 648 ico_l = coset(lxa, lya, lza) 649 jco_l = coset(MAX(lxb - 1, 0), (lyb + 1), lzb) 650 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 651 652 ico_l = coset(lxa, lya, lza) 653 jco_l = coset((lxb + 1), (lyb + 1), lzb) 654 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 655 656 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 657 jco_l = coset(lxb, (lyb + 1), lzb) 658 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lxa*pab(o1 + ico, o2 + jco) 659 660 ico_l = coset((lxa + 1), lya, lza) 661 jco_l = coset(lxb, (lyb + 1), lzb) 662 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 663 664 ELSEIF (idir == 1 .AND. ir == 3) THEN 665 666 ico_l = coset(lxa, lya, lza) 667 jco_l = coset(MAX(lxb - 1, 0), lyb, (lzb + 1)) 668 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 669 670 ico_l = coset(lxa, lya, lza) 671 jco_l = coset((lxb + 1), lyb, (lzb + 1)) 672 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 673 674 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 675 jco_l = coset(lxb, lyb, (lzb + 1)) 676 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lxa*pab(o1 + ico, o2 + jco) 677 678 ico_l = coset((lxa + 1), lya, lza) 679 jco_l = coset(lxb, lyb, (lzb + 1)) 680 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 681 682 ELSEIF (idir == 2 .AND. ir == 1) THEN 683 684 ico_l = coset(lxa, lya, lza) 685 jco_l = coset((lxb + 1), MAX(lyb - 1, 0), lzb) 686 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 687 688 ico_l = coset(lxa, lya, lza) 689 jco_l = coset((lxb + 1), (lyb + 1), lzb) 690 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 691 692 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 693 jco_l = coset((lxb + 1), lyb, lzb) 694 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lya*pab(o1 + ico, o2 + jco) 695 696 ico_l = coset(lxa, (lya + 1), lza) 697 jco_l = coset((lxb + 1), lyb, lzb) 698 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 699 700 ELSEIF (idir == 2 .AND. ir == 2) THEN 701 702 ico_l = coset(lxa, lya, lza) 703 jco_l = coset(lxb, lyb, lzb) 704 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 705 706 ico_l = coset(lxa, lya, lza) 707 jco_l = coset(lxb, (lyb + 2), lzb) 708 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 709 710 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 711 jco_l = coset(lxb, (lyb + 1), lzb) 712 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lya*pab(o1 + ico, o2 + jco) 713 714 ico_l = coset(lxa, (lya + 1), lza) 715 jco_l = coset(lxb, (lyb + 1), lzb) 716 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 717 718 ELSEIF (idir == 2 .AND. ir == 3) THEN 719 720 ico_l = coset(lxa, lya, lza) 721 jco_l = coset(lxb, MAX(lyb - 1, 0), (lzb + 1)) 722 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 723 724 ico_l = coset(lxa, lya, lza) 725 jco_l = coset(lxb, (lyb + 1), (lzb + 1)) 726 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 727 728 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 729 jco_l = coset(lxb, lyb, (lzb + 1)) 730 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lya*pab(o1 + ico, o2 + jco) 731 732 ico_l = coset(lxa, (lya + 1), lza) 733 jco_l = coset(lxb, lyb, (lzb + 1)) 734 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 735 736 ELSEIF (idir == 3 .AND. ir == 1) THEN 737 738 ico_l = coset(lxa, lya, lza) 739 jco_l = coset((lxb + 1), lyb, MAX(lzb - 1, 0)) 740 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 741 742 ico_l = coset(lxa, lya, lza) 743 jco_l = coset((lxb + 1), lyb, (lzb + 1)) 744 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 745 746 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 747 jco_l = coset((lxb + 1), lyb, lzb) 748 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lza*pab(o1 + ico, o2 + jco) 749 750 ico_l = coset(lxa, lya, (lza + 1)) 751 jco_l = coset((lxb + 1), lyb, lzb) 752 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 753 754 ELSEIF (idir == 3 .AND. ir == 2) THEN 755 756 ico_l = coset(lxa, lya, lza) 757 jco_l = coset(lxb, (lyb + 1), MAX(lzb - 1, 0)) 758 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 759 760 ico_l = coset(lxa, lya, lza) 761 jco_l = coset(lxb, (lyb + 1), (lzb + 1)) 762 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 763 764 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 765 jco_l = coset(lxb, (lyb + 1), lzb) 766 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lza*pab(o1 + ico, o2 + jco) 767 768 ico_l = coset(lxa, lya, (lza + 1)) 769 jco_l = coset(lxb, (lyb + 1), lzb) 770 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 771 772 ELSEIF (idir == 3 .AND. ir == 3) THEN 773 774 ico_l = coset(lxa, lya, lza) 775 jco_l = coset(lxb, lyb, lzb) 776 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 777 778 ico_l = coset(lxa, lya, lza) 779 jco_l = coset(lxb, lyb, (lzb + 2)) 780 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 781 782 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 783 jco_l = coset(lxb, lyb, (lzb + 1)) 784 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lza*pab(o1 + ico, o2 + jco) 785 786 ico_l = coset(lxa, lya, (lza + 1)) 787 jco_l = coset(lxb, lyb, (lzb + 1)) 788 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 789 790 END IF 791 792 END SUBROUTINE prepare_ardb_m_darb 793 794END MODULE grid_modify_pab_block 795