1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18c /* deck cclr_fa */ 19*=====================================================================* 20 SUBROUTINE CCLR_FA ( LABELA, ISYMTA, ! inp: label/symmetry A 21 & LISTB, ITAMPB, ! inp: B resp. amplit. 22 & LISTC, IZETVC, ! inp: C resp. zeta vec. 23 & XINT, ! inp: integrals of operator 24 & WORK, LWORK )! work space 25*---------------------------------------------------------------------* 26* 27* Purpose: transformation of a response vector with a F matrix 28* where the hamiltonian has been substituted by a 29* perturbation operator 30* 31* F^C{A} * t^B = <lambda^C|[[A,t^B],tau_nu]|CC> 32* 33* JK+OC. CCMM: Allow for input of integrals if 34* LABELA .eq. 'GIVE INT' 35* 36* symmetries/variables: 37* 38* ISYRES : result vector GAMMA1, GAMMA2 39* ISYCTR : lagrangian multipliers (zeta vector) CTR1, CTR2 40* ISYMTA : A perturbation 41* ISYMTB : B response amplitudes 42* 43* Note: the single and double excitation parts of the result GAMMA2 44* are returned at the beginning of the work space in 45* WORK(1)... WORK(NT1AM(ISYRES)) 46* WORK(NT1AM(ISYRES)+1)... WORK(NT1AM(ISYRES)+NT2AM(ISYRES)) 47* (double excitation part will be stored in packed form) 48* 49* Written by Christof Haettig, October 1996. 50* 51*=====================================================================* 52#if defined (IMPLICIT_NONE) 53 IMPLICIT NONE 54#else 55# include "implicit.h" 56#endif 57#include "priunit.h" 58#include "iratdef.h" 59#include "cbieri.h" 60#include "mxcent.h" 61#include "eribuf.h" 62#include "maxorb.h" 63#include "distcl.h" 64#include "ccorb.h" 65#include "ccisao.h" 66#include "ccsdsym.h" 67#include "ccsdinp.h" 68 69* local parameters: 70 CHARACTER MSGDBG*(17) 71 PARAMETER (MSGDBG='[debug] CCLR_FA> ') 72#if defined (SYS_CRAY) 73 REAL ONE, TWO, THREE 74#else 75 DOUBLE PRECISION ONE, TWO, THREE 76#endif 77 PARAMETER (ONE = 1.0d0, TWO = 2.0d0, THREE = 3.0d0) 78 LOGICAL LOCDBG 79 PARAMETER (LOCDBG = .FALSE.) 80 INTEGER KDUM 81 PARAMETER (KDUM = +99 999 999) ! dummy address 82 83 CHARACTER*8 LABELA 84 CHARACTER*10 MODEL 85 CHARACTER LISTB*(*), LISTC*(*) 86 INTEGER ISYRES, ISYCTR, ISYMTA, ISYMTB, LWORK 87 INTEGER ITAMPB, IZETVC 88 89 INTEGER ISYTATB, ISYMJ, ISYMB, ISYMXY, ISYMI, ISYMA, IRREP, ISYM 90 INTEGER KBTAOO, KBTAVV, KPERTA, KT1AMPB, KT2AMPB, KCTR1, KCTR2 91 INTEGER KCTMO, KT1AMP0, KLAMDP0, KLAMDH0, KOFF1, KOFF2, KSCR 92 INTEGER KEND1, KEND2, KEND3, KEND1A, KXBMAT, KYBMAT, IERR 93 INTEGER LEND1, LEND2, LEND3, LEND1A, KGAMMA1, KGAMMA2, KEND0 94 INTEGER IOPT, MAXJ, NIJ, NJI, NAB, NBA, KEMAT1, KEMAT2 95 96#if defined (SYS_CRAY) 97 REAL FREQB 98 REAL SWAP 99 REAL WORK(LWORK), XINT(*) 100#else 101 DOUBLE PRECISION FREQB 102 DOUBLE PRECISION SWAP 103 DOUBLE PRECISION WORK(LWORK), XINT(*) 104#endif 105 106 INTEGER ILSTSYM 107 108*---------------------------------------------------------------------* 109* begin: 110*---------------------------------------------------------------------* 111 IF (CCSDT) THEN 112 WRITE (LUPRI,'(/1x,a)') 'F{A} matrix transformations not ' 113 & //'implemented for triples yet...' 114 CALL QUIT('Triples not implemented for G '// 115 & 'matrix transformations') 116 END IF 117 118 IF ( .not. (CCS .or. CC2 .or. CCSD) ) THEN 119 WRITE (LUPRI,'(/1x,a)') 'CCLR_FA called for a Coupled Cluster ' 120 & //'method not implemented in CCLR_FA...' 121 CALL QUIT('Unknown CC method in CCLR_FA.') 122 END IF 123 124*---------------------------------------------------------------------* 125* set & check symmetries: 126*---------------------------------------------------------------------* 127 ISYMTB = ILSTSYM(LISTB,ITAMPB) ! B 128 ISYCTR = ILSTSYM(LISTC,IZETVC) ! C 129 ISYTATB = MULD2H(ISYMTA,ISYMTB) ! A x B 130 ISYRES = MULD2H(ISYCTR,ISYTATB) ! A x B x C 131 132 IF (LOCDBG) THEN 133 WRITE (LUPRI,*) 'LISTB,ITAMPB,ISYMTB:',LISTB,ITAMPB,ISYMTB 134 WRITE (LUPRI,*) 'LISTC,IZETVC,ISYCTR:',LISTC,IZETVC,ISYCTR 135 END IF 136 137 138 IF (ISYMOP .NE. 1) THEN 139 WRITE (LUPRI,'(/1x,a)') 'non-total-symmetric MO integrals?!... ' 140 & //'CCLR_G has never been debugged for that!...' 141 END IF 142 143 IF (MULD2H(ISYCTR,ISYTATB) .NE. ISYRES) THEN 144 CALL QUIT('Symmetry mismatch in CCLR_FA.') 145 END IF 146 147*---------------------------------------------------------------------* 148* flush print unit 149*---------------------------------------------------------------------* 150 Call FLSHFO(LUPRI) 151 152 IF (LOCDBG) THEN 153 WRITE (LUPRI,'(/1x,a,i15)') 'work space in CCLR_FA:',LWORK 154 END IF 155*---------------------------------------------------------------------* 156* initialize pointer for work space and allocate memory for 157* 1) single excitation part of the result vector 158* 2) one-index transformed perturbation integrals A^B (occ/occ block) 159* 3) one-index transformed perturbation integrals A^B (vir/vir block) 160* 4) perturbation integrals A 161* 5) singles part of response amplitudes T1^B 162* 6) singles part of zeroth order lagrangian multipliers 163*---------------------------------------------------------------------* 164 KGAMMA1 = 1 165 KBTAOO = KGAMMA1 + NT1AM(ISYRES) 166 KBTAVV = KBTAOO + NMATIJ(ISYTATB) 167 KPERTA = KBTAVV + NMATAB(ISYTATB) 168 KT1AMPB = KPERTA + NT1AM(ISYMTA) 169 KCTR1 = KT1AMPB + NT1AM(ISYMTB) 170 KEND1 = KCTR1 + NT1AM(ISYCTR) 171 LEND1 = LWORK - KEND1 172 173 IF (LEND1 .LT. 0) THEN 174 CALL QUIT('Insufficient work space in CCLR_FA.') 175 END IF 176 177*---------------------------------------------------------------------* 178* initialize single excitation part of result vector GAMMA1: 179*---------------------------------------------------------------------* 180 Call DZERO (WORK(KGAMMA1), NT1AM(ISYRES)) 181 182*---------------------------------------------------------------------* 183* for CCS and zeroth-order zeta vector all contributions vanish: 184*---------------------------------------------------------------------* 185 IF (CCS .AND. LISTC(1:2).EQ.'L0') RETURN 186 187*---------------------------------------------------------------------* 188* read singles parts for B response amplitudes and zeta vector: 189*---------------------------------------------------------------------* 190 IOPT = 1 191 CALL CC_RDRSP(LISTB,ITAMPB,ISYMTB,IOPT,MODEL, 192 & WORK(KT1AMPB),WORK(KDUM) ) 193 194 195 IOPT = 1 196 Call CC_RDRSP(LISTC,IZETVC,ISYCTR,IOPT,MODEL, 197 & WORK(KCTR1),WORK(KDUM)) 198 199 IF (LOCDBG) THEN 200 CAll AROUND('response T amplitudes B:') 201 WRITE (LUPRI,*) 'LIST/INDEX:',LISTB,ITAMPB 202 WRITE (LUPRI,*) 'Symmetry: ',ISYMTB 203 CAll CC_PRP(WORK(KT1AMPB),WORK(KDUM),ISYMTB,1,0) 204 CALL AROUND('CC lagrange multipliers') 205 CALL CC_PRP(WORK(KCTR1), WORK(KDUM), ISYCTR, 1, 0) 206 END IF 207 208*---------------------------------------------------------------------* 209* read & resort one-electron integrals for operator A: 210*---------------------------------------------------------------------* 211 KCTMO = KEND1 212 KT1AMP0 = KCTMO + N2BST(ISYMTA) 213 KLAMDP0 = KT1AMP0 + NT1AM(ISYMOP) 214 KLAMDH0 = KLAMDP0 + NLAMDT 215 KEND1A = KLAMDH0 + NLAMDT 216 LEND1A = LWORK - KEND1A 217 218 IF (LEND1A .LT. 0) THEN 219 CALL QUIT('Insufficient work space in CCLR_FA.') 220 END IF 221C 222C JK+OC, CCSLV 223 IF (LABELA.EQ.'GIVE INT') THEN 224 CALL DCOPY(N2BST(ISYMTA),XINT(1),1,WORK(KCTMO),1) 225 ELSE 226* read the AO integrals: 227 CALL CCPRPAO(LABELA,.TRUE.,WORK(KCTMO),IRREP,ISYM,IERR, 228 & WORK(KEND1A),LEND1A) 229 IF (IERR.NE.0 .OR. IRREP.NE.ISYMTA) THEN 230 CALL QUIT('CCLR_FA: error while reading operator '//LABELA) 231 END IF 232 END IF 233C 234* get MO coefficients: 235 CALL DZERO(WORK(KT1AMP0),NT1AMX) 236 CALL LAMMAT(WORK(KLAMDP0),WORK(KLAMDH0),WORK(KT1AMP0), 237 & WORK(KEND1A),LEND1A) 238 239* transform one-electron integrals in place: 240 CALL CC_FCKMO(WORK(KCTMO),WORK(KLAMDP0),WORK(KLAMDH0), 241 & WORK(KEND1A),LEND1A,ISYMTA,1,1) 242 243* resort occupied/virtual block to T1 like storage: 244 CALL DZERO(WORK(KPERTA),NT1AM(ISYMTA)) 245 DO ISYMJ = 1, NSYM 246 ISYMB = MULD2H(ISYMJ,ISYMTA) 247 248 DO J = 1, NRHF(ISYMJ) 249 DO B = 1, NVIR(ISYMB) 250 KOFF1 = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1) + B 251 KOFF2 = IFCVIR(ISYMJ,ISYMB) + NORB(ISYMJ)*(B-1) + J 252 253 WORK(KPERTA-1+KOFF1) = WORK(KCTMO-1+KOFF2) 254 END DO 255 END DO 256 END DO 257 258 IF (LOCDBG) THEN 259 WRITE (LUPRI,*) MSGDBG,' A integrals in MO basis:' 260 WRITE (LUPRI,*) MSGDBG,' label, symmetry:',LABELA,ISYMTA 261 Call CC_PRP(WORK(KPERTA),WORK(KDUM),ISYMTA,1,0) 262 END IF 263 264*---------------------------------------------------------------------* 265* calculate A perturbation integrals one-index transformed with 266* the B response amplitudes T1^B: 267*---------------------------------------------------------------------* 268* occ/occ block: 269 Call CCG_1ITROO(WORK(KBTAOO), ISYTATB, 270 & WORK(KPERTA), ISYMTA, 271 & WORK(KT1AMPB),ISYMTB ) 272 273* vir/vir block: 274 Call CCG_1ITRVV(WORK(KBTAVV), ISYTATB, 275 & WORK(KPERTA), ISYMTA, 276 & WORK(KT1AMPB),ISYMTB ) 277 278*=====================================================================* 279* CCS part: < Zeta_1 | [tA^B, tau_1] | HF> 280*=====================================================================* 281* do one-index transformation with Zeta vector: 282 IOPT = 2 283 Call CCG_1ITRVO(WORK(KGAMMA1),ISYRES,WORK(KBTAOO),WORK(KBTAVV), 284 & ISYTATB,WORK(KCTR1),ISYCTR,IOPT ) 285 286 IF (LOCDBG) THEN 287 WRITE (LUPRI,*) MSGDBG,'one-index trans. A (occ/occ block):' 288 WRITE (LUPRI,'(5f12.6)') (WORK(KBTAOO-1+I),I=1,NMATIJ(ISYTATB)) 289 WRITE (LUPRI,*) MSGDBG,'one-index trans. A (vir/vir block):' 290 WRITE (LUPRI,'(2f12.6)') (WORK(KBTAVV-1+I),I=1,NMATAB(ISYTATB)) 291 WRITE (LUPRI,*) MSGDBG, 292 * 'contrib. of one-index trans. A to GAMMA:' 293 Call CC_PRP(WORK(KGAMMA1),WORK(KDUM),ISYRES,1,0) 294 END IF 295 296*---------------------------------------------------------------------* 297* end of CCS part 298*---------------------------------------------------------------------* 299 300 If (CCS) Return 301 302*=====================================================================* 303* CC2/CCSD part for the singles: <Zeta_2| [[A,T2B], tau_1] |HF> 304*=====================================================================* 305 306*---------------------------------------------------------------------* 307* memory allocation: 308* 1) double excitation part of response amplitudes T2B (packed) 309* 2) double excitation part of zeta vector (squared) 310* 3) double excitation part of zeta vector (packed) 311* N.B. we account here for the fact, that the packed double excitation 312* part of the result vector will be returned at the beginning of the 313* work space, so we make sure, that there is enough space before 314* the zeta vector to store there later on GAMMA2 315*---------------------------------------------------------------------* 316 KT2AMPB = KEND1 317 KCTR2 = KT2AMPB + MAX( NT2AM(ISYMTB), NT2AM(ISYRES) ) 318 KEND2 = KCTR2 + NT2SQ(ISYCTR) 319 LEND2 = LWORK - KEND2 320 321 IF (LEND2 .LT. NT2AM(ISYCTR) ) THEN 322 CALL QUIT('Insufficient work space in CCLR_FA.') 323 END IF 324 325*---------------------------------------------------------------------* 326* read response amplitudes T2B and scale the diagonal: 327*---------------------------------------------------------------------* 328 IOPT = 2 329 CALL CC_RDRSP(LISTB,ITAMPB,ISYMTB,IOPT,MODEL, 330 & WORK(KDUM),WORK(KT2AMPB) ) 331 332 CAll CCLR_DIASCL(WORK(KT2AMPB),TWO,ISYMTB) 333 334 IF (LOCDBG) THEN 335 WRITE (LUPRI,*) MSGDBG, 'B response amplitudes:' 336 Call CC_PRP(WORK(KT1AMPB),WORK(KT2AMPB),ISYMTB,1,1) 337 END IF 338 339*---------------------------------------------------------------------* 340* read packed lagrangian multipliers and square them up: 341*---------------------------------------------------------------------* 342 IOPT = 2 343 Call CC_RDRSP(LISTC,IZETVC,ISYCTR,IOPT,MODEL, 344 & WORK(KDUM),WORK(KEND2)) 345 346 CALL CC_T2SQ (WORK(KEND2), WORK(KCTR2), ISYCTR) 347 348*---------------------------------------------------------------------* 349* calculate X^B and Y^B intermediates: 350*---------------------------------------------------------------------* 351 ISYMXY = MULD2H(ISYCTR,ISYMTB) 352 353 KXBMAT = KEND2 354 KYBMAT = KXBMAT + NMATIJ(ISYMXY) 355 KSCR = KYBMAT + NMATAB(ISYMXY) 356 KEND3 = KSCR + NT1AM(ISYRES) 357 LEND3 = LWORK - KEND3 358 359 If (LEND3 .LT. 0) THEN 360 CALL QUIT('Insufficient work space in CCLR_FA.') 361 END IF 362 363* calculate X^C & Y^C intermediate: 364 Call CC_XI(WORK(KXBMAT),WORK(KCTR2), ISYCTR, 365 & WORK(KT2AMPB),ISYMTB,WORK(KEND3),LEND3) 366 367 Call CC_YI(WORK(KYBMAT),WORK(KCTR2), ISYCTR, 368 & WORK(KT2AMPB),ISYMTB,WORK(KEND3),LEND3) 369 370 IF (LOCDBG) THEN 371 WRITE (LUPRI,*) MSGDBG,'response X intermediate:' 372 WRITE (LUPRI,'(5f12.6)') (WORK(KXBMAT-1+I),I=1,NMATIJ(ISYMXY)) 373 WRITE (LUPRI,*) MSGDBG,'response Y intermediate:' 374 WRITE (LUPRI,'(2f12.6)') (WORK(KYBMAT-1+I),I=1,NMATAB(ISYMXY)) 375 END IF 376 377* calculate XY^B: XY^B_ij = X^B_ji, XY^B_bd = -Y^B_bd 378* 1.) transpose X^B intermediate 379 DO ISYMI = 1, NSYM 380 ISYMJ = MULD2H(ISYMI,ISYMXY) 381 IF (ISYMJ .LE. ISYMI) THEN 382 DO I = 1, NRHF(ISYMI) 383 MAXJ = NRHF(ISYMJ) 384 IF (ISYMJ .EQ. ISYMI) MAXJ = I-1 385 DO J = 1, MAXJ 386 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 387 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I-1) + J 388 SWAP = WORK(KXBMAT-1+NIJ) 389 WORK(KXBMAT-1+NIJ) = WORK(KXBMAT-1+NJI) 390 WORK(KXBMAT-1+NJI) = SWAP 391 END DO 392 END DO 393 END IF 394 END DO 395 396* 2.) multiply Y^B intermediate with -1: 397 Call DSCAL(NMATAB(ISYMXY), -ONE, WORK(KYBMAT), 1) 398 399 400* do one-index transformation of XY^B with A integrals: 401 IOPT = 2 402 Call CCG_1ITRVO(WORK(KSCR),ISYRES, 403 & WORK(KXBMAT),WORK(KYBMAT),ISYMXY, 404 & WORK(KPERTA),ISYMTA, IOPT ) 405 406* add contribution to GAMMA1: 407 Call DAXPY (NT1AM(ISYRES), ONE, WORK(KSCR),1, WORK(KGAMMA1), 1) 408 409 IF (LOCDBG) THEN 410 WRITE (LUPRI,*) MSGDBG,'A integrals:' 411 WRITE (LUPRI,'(5f12.6)') (WORK(KPERTA-1+I),I=1,NT1AM(ISYMTA)) 412 WRITE (LUPRI,*) MSGDBG, 'CC2/CCSD contribution to singles part:' 413 Call CC_PRP(WORK(KSCR),WORK(KDUM),ISYRES,1,0) 414 WRITE (LUPRI,*) MSGDBG, 'GAMMA1 now:' 415 Call CC_PRP(WORK(KGAMMA1),WORK(KDUM),ISYRES,1,0) 416 END IF 417 418*---------------------------------------------------------------------* 419 420*=====================================================================* 421* CC2/CCSD part for the doubles: <Zeta_2| [[A,T1B], tau_2] |HF> 422*=====================================================================* 423 424*---------------------------------------------------------------------* 425* reorganize work space, so that the result vector GAMMA2 can be 426* stored at the early beginning of the work space 427*---------------------------------------------------------------------* 428 KGAMMA1 = 1 429 KGAMMA2 = KGAMMA1 + NT1AM(ISYRES) 430 KEND0 = KGAMMA2 + NT2AM(ISYRES) 431 432 IF (KEND0 .GT. KCTR2) THEN 433 CALL QUIT('memory organization mixed up in CCLR_FA.') 434 END IF 435 436 KEMAT1 = KCTR2 + NT2SQ(ISYCTR) 437 KEMAT2 = KEMAT1 + NMATAB(ISYTATB) 438 KEND3 = KEMAT2 + NMATIJ(ISYTATB) 439 LEND3 = LWORK - KEND3 440 441 IF ( LEND3 .LT. 0 ) THEN 442 CALL QUIT('Insufficient work space in CCLR_FA.') 443 END IF 444 445*---------------------------------------------------------------------* 446* transpose tA^B(a b) --> EMAT1(b a) 447*---------------------------------------------------------------------* 448 DO ISYMA = 1, NSYM 449 ISYMB = MULD2H(ISYMA,ISYTATB) 450 DO A = 1, NVIR(ISYMA) 451 DO B = 1, NVIR(ISYMB) 452 NAB = IMATAB(ISYMA,ISYMB) + NVIR(ISYMA)*(B-1) + A 453 NBA = IMATAB(ISYMB,ISYMA) + NVIR(ISYMB)*(A-1) + B 454 455 WORK(KEMAT1 - 1 + NBA) = WORK(KBTAVV - 1 + NAB) 456 END DO 457 END DO 458 END DO 459 460 461*---------------------------------------------------------------------* 462* transpose tA^B(i j) --> EMAT2(j i) 463*---------------------------------------------------------------------* 464 DO ISYMI = 1, NSYM 465 ISYMJ = MULD2H(ISYMI,ISYTATB) 466 DO I = 1, NRHF(ISYMI) 467 DO J = 1, NRHF(ISYMJ) 468 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J-1) + I 469 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I-1) + J 470 471 WORK(KEMAT2 - 1 + NJI) = WORK(KBTAOO - 1 + NIJ) 472 END DO 473 END DO 474 END DO 475 476 IF (LOCDBG) THEN 477 WRITE (LUPRI,*) MSGDBG,'one-index trans. A (occ/occ block):' 478 WRITE (LUPRI,'(5f12.6)') (WORK(KBTAOO-1+I),I=1,NMATIJ(ISYTATB)) 479 WRITE (LUPRI,*) MSGDBG,'one-index trans. A (vir/vir block):' 480 WRITE (LUPRI,'(2f12.6)') (WORK(KBTAVV-1+I),I=1,NMATAB(ISYTATB)) 481 WRITE (LUPRI,*) MSGDBG,'EMAT2:' 482 WRITE (LUPRI,'(5f12.6)') (WORK(KEMAT2-1+I),I=1,NMATIJ(ISYTATB)) 483 WRITE (LUPRI,*) MSGDBG,'EMAT1:' 484 WRITE (LUPRI,'(2f12.6)') (WORK(KEMAT1-1+I),I=1,NMATAB(ISYTATB)) 485 END IF 486 487*---------------------------------------------------------------------* 488* combine EMAT1/EMAT2 with lagrangian multipliers: 489* (note: this overwrites the intermedites stored at the beginning 490* of the work space...) 491*---------------------------------------------------------------------* 492* initialize GAMMA2: 493 CALL DZERO(WORK(KGAMMA2),NT2AM(ISYRES)) 494 495* do the caculation: 496 Call CCRHS_E(WORK(KGAMMA2),WORK(KCTR2),WORK(KEMAT1), 497 & WORK(KEMAT2), WORK(KEND3), LEND3, ISYCTR, ISYTATB) 498 499*---------------------------------------------------------------------* 500 IF (LOCDBG) THEN 501 WRITE (LUPRI,*) MSGDBG, 'GAMMA:' 502 Call CC_PRP(WORK(KGAMMA1),WORK(KGAMMA2),ISYRES,1,1) 503 END IF 504C 505 RETURN 506 END 507*=====================================================================* 508* END OF SUBROUTINE CCLR_FA * 509*=====================================================================* 510