1 subroutine sp_start(lout,lfntop,filtop,lfnrst,filrst, 2 + lsyn,filsyn,nsyn, 3 + rsht,rlng,rinp,rsgmi, 4 + npxi,npyi,npzi,nbxi,nbyi,nbzi, 5 + npbt,nbxt,boxt,vlatt,lpbc, 6 + nwmi,mwmi,nwai,mwai,nsfi,msfi,nsmi,msmi,nsai,msai, 7 + ldb,lbp,fld,lpol,lfre,temp,tempw,temps,lqd,iguidi, 8 + ldbg,idbg,prjct,mbbi,nseri,isld,irset,ictrl,nseqi, 9 + i_lseqi,ndumsi,nbgeti,npreci,madbxi) 10c 11c 12c $Id$ 13c 14 implicit none 15c 16#include "sp_common.fh" 17#include "mafdecls.fh" 18#include "msgids.fh" 19#include "global.fh" 20c 21 logical sp_rdrest,sp_diffbb 22 external sp_rdrest,sp_diffbb 23c 24 integer lfntop,lfnrst,lsyn,nsyn,ldbg,idbg,lout 25 character*255 filtop,filrst,filsyn 26 real*8 rsht,rlng,rinp,boxt(3),vlatt(3,3),rtemp(3) 27 real*8 temp,tempw,temps,fld,rsgmi 28 integer itemp(2) 29 integer npxi,npyi,npzi,nbxi,nbyi,nbzi,isld,irset,nseqi,i_lseqi 30 integer nwmi,mwmi,nwai,mwai,nsfi,msfi,nsmi,msmi,nsai,msai 31 integer lenscr,ldb,lbp,npbt,nbxt,iguidi,mbbi,nseri,ictrl,ndumsi 32 integer nbgeti,npreci,madbxi 33 logical lpol,lfre,ignore,lpbc,lqd 34 character*80 prjct 35c 36 integer i,j,l_f,i_f 37c 38 if(.not.ma_verify_allocator_stuff()) then 39 call md_abort('ERROR IN MEMORY USE AT START SP_START',0) 40 endif 41c 42 me=ga_nodeid() 43 np=ga_nnodes() 44c 45 project=prjct 46c 47 idebug=idbg 48 lfndbg=ldbg 49 icntrl=ictrl 50 lfnout=lout 51c 52 lfnsyn=lsyn 53 nfsync=nsyn 54 isload=isld 55 ireset=irset 56 nbget=nbgeti 57 nprec=npreci 58 ibget=nbget 59 madbox=madbxi 60c 61 rshort=rsht 62 rlong=max(rsht,rlng) 63 rbox=rinp 64 ntwin=1 65 if(rlng.gt.rsht) ntwin=2 66 ltwin=rlng.gt.rsht 67 lnode0=lqmd 68 lqmd=lqd 69c 70 loadb=ldb 71 lbpair=lbp 72 factld=fld 73c 74 lpola=lpol 75 lfree=lfre 76c 77 nbx=nbxi 78 nby=nbyi 79 nbz=nbzi 80 nbxin=nbxi 81 nbyin=nbyi 82 nbzin=nbzi 83c 84 npx=npxi 85 npy=npyi 86 npz=npzi 87c 88 mbbl=0 89 nbbdif=-1 90c 91 mwm=mwmi 92 msa=max(msai,msmi) 93 mbbreq=mbbi 94 mbblp=0 95 nserie=nseri 96c 97 iguide=iguidi 98c 99 npack=0 100 npackw=0 101c 102 rsgm=rsgmi 103c 104 call sp_nrnode() 105c 106 call sp_dimens(lfnrst,filrst) 107c 108 call sp_alloc() 109c 110 call sp_decomp(int_mb(i_iown),dbl_mb(i_boxs), 111 + int_mb(i_buren),int_mb(i_bindex)) 112c 113 call sp_initip(int_mb(i_iown),int_mb(i_ipl)) 114c 115 call sp_numbb(int_mb(i_iown),dbl_mb(i_boxs)) 116c 117 if(sp_diffbb(dbl_mb(i_boxs),int_mb(i_rng))) 118 + call sp_listbb(int_mb(i_iown),dbl_mb(i_boxs),int_mb(i_bb)) 119c 120 call sp_alloc2() 121c 122 if(ireset.eq.0) then 123 ignore=sp_rdrest(lfnrst,filrst,dbl_mb(i_boxs)) 124 endif 125 if(np.gt.1) then 126 call ga_brdcst(msp_02,nsm,ma_sizeof(mt_int,1,mt_byte),0) 127 endif 128 msm=max(1,nsm) 129 if(.not.ma_push_get(mt_dbl,msm*3,'xscr',l_xscr,i_xscr)) 130 + call md_abort('Failed to allocate xscr',0) 131c 132 lenscr=3*max(mwm*mwa,msa) 133 if(.not.ma_push_get(mt_dbl,lenscr,'x',l_x,i_x)) 134 + call md_abort('Failed to allocate x',0) 135 if(.not.ma_push_get(mt_dbl,lenscr,'v',l_v,i_v)) 136 + call md_abort('Failed to allocate v',0) 137 if(.not.ma_push_get(mt_dbl,lenscr,'f',l_f,i_f)) 138 + call md_abort('Failed to allocate f',0) 139 if(.not.ma_push_get(mt_dbl,3*mwm*mwa,'r',l_r,i_r)) 140 + call md_abort('Failed to allocate r',0) 141 lenscr=max(miw2*mwm,mis2*msa) 142 if(.not.ma_push_get(mt_int,lenscr,'i',l_i,i_i)) 143 + call md_abort('Failed to allocate i',0) 144 lenscr=ma_inquire_avail(mt_byte)/ 145 + ((9*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+ 146 + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1 147 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx)) 148 + call md_abort('Failed to allocate bx',0) 149 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv)) 150 + call md_abort('Failed to allocate bv',0) 151 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bf',l_bf,i_bf)) 152 + call md_abort('Failed to allocate bf',0) 153 if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br)) 154 + call md_abort('Failed to allocate br',0) 155 if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi)) 156 + call md_abort('Failed to allocate bi',0) 157 if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n)) 158 + call md_abort('Failed to allocate n',0) 159c 160 call ga_sync() 161c 162 if(.not.ma_verify_allocator_stuff()) then 163 call md_abort('ERROR IN MEMORY USE BEFORE RDRST IN SP_START',0) 164 endif 165c 166 call sp_rdrst(lfnrst,filrst,lfntop,filtop, 167 + temp,tempw,temps,int_mb(i_ipl), 168 + dbl_mb(i_x),dbl_mb(i_v),dbl_mb(i_f),dbl_mb(i_r),int_mb(i_i), 169 + dbl_mb(i_x),dbl_mb(i_v),dbl_mb(i_f),dbl_mb(i_xscr),int_mb(i_i), 170 + dbl_mb(i_bx),dbl_mb(i_bv),dbl_mb(i_bf),dbl_mb(i_br), 171 + int_mb(i_bi),lenscr, 172 + dbl_mb(i_bx),dbl_mb(i_bv),dbl_mb(i_bf),int_mb(i_bi), 173 + lenscr,int_mb(i_n), 174 + int_mb(i_iown),dbl_mb(i_boxs),int_mb(i_lseq),int_mb(i_sndx)) 175c 176 if(.not.ma_verify_allocator_stuff()) then 177 call md_abort('ERROR IN MEMORY USE AFTER RDRST IN SP_START',0) 178 endif 179c 180 if(.not.ma_pop_stack(l_n)) 181 + call md_abort('Failed to deallocate n',0) 182 if(.not.ma_pop_stack(l_bi)) 183 + call md_abort('Failed to deallocate bi',0) 184 if(.not.ma_pop_stack(l_br)) 185 + call md_abort('Failed to deallocate br',0) 186 if(.not.ma_pop_stack(l_bf)) 187 + call md_abort('Failed to deallocate bf',0) 188 if(.not.ma_pop_stack(l_bv)) 189 + call md_abort('Failed to deallocate bv',0) 190 if(.not.ma_pop_stack(l_bx)) 191 + call md_abort('Failed to deallocate bx',0) 192 if(.not.ma_pop_stack(l_i)) 193 + call md_abort('Failed to deallocate i',0) 194 if(.not.ma_pop_stack(l_r)) 195 + call md_abort('Failed to deallocate r',0) 196 if(.not.ma_pop_stack(l_f)) 197 + call md_abort('Failed to deallocate f',0) 198 if(.not.ma_pop_stack(l_v)) 199 + call md_abort('Failed to deallocate v',0) 200 if(.not.ma_pop_stack(l_x)) 201 + call md_abort('Failed to deallocate x',0) 202c 203 if(np.gt.1) then 204 rtemp(1)=temp 205 rtemp(2)=tempw 206 rtemp(3)=temps 207 call ga_brdcst(msp_01,rtemp,ma_sizeof(mt_dbl,3,mt_byte),0) 208 temp=rtemp(1) 209 tempw=rtemp(2) 210 temps=rtemp(3) 211 itemp(1)=nsm 212 itemp(2)=nsf 213 call ga_brdcst(msp_02,itemp,ma_sizeof(mt_int,2,mt_byte),0) 214 nsm=itemp(1) 215 nsf=itemp(2) 216 endif 217c 218 nwmi=nwm 219 nwai=nwa 220 nsmi=nsm 221 nsai=nsa 222c 223 mwmi=mwm 224 mwai=mwa 225 msai=msa 226c 227 npbt=npbtyp 228 nbxt=nbxtyp 229c 230 do 1 i=1,3 231 boxt(i)=box(i) 232 do 2 j=1,3 233 vlatt(i,j)=vlat(i,j) 234 2 continue 235 1 continue 236c 237 nsfi=nsf 238 msfi=max(1,nsf) 239 msmi=max(1,nsm) 240c 241 third=1.0d0/3.0d0 242 nldup=-2 243 ipairf=-1 244 ipairt=-1 245 lpipo=.false. 246c 247 lpbc=npbtyp.gt.0 248c 249 if(me.eq.0) then 250 if(nfsync.gt.0) then 251 open(unit=lfnsyn,file=filsyn(1:index(filsyn,' ')-1), 252 + status='unknown') 253 write(lfnsyn,3000) np 254 3000 format(i5) 255 endif 256 endif 257c 258 nseqi=nseq 259 i_lseqi=i_lseq 260c 261 ndumsi=ndums 262c 263 if(.not.ma_verify_allocator_stuff()) then 264 call md_abort('ERROR IN MEMORY USE AT EXIT SP_START',0) 265 endif 266c 267 return 268 end 269 subroutine sp_dimens(lfnrst,filrst) 270c 271 implicit none 272c 273#include "sp_common.fh" 274#include "msgids.fh" 275#include "mafdecls.fh" 276#include "global.fh" 277#include "geom.fh" 278#include "util.fh" 279c 280 integer lfnrst 281 character*255 filrst 282c 283 character*1 cdum 284 integer i,j,jdum,ibx,iby,ibz,itemp(11) 285 real*8 rtemp(5),rsgmr 286c 287 mbbl=0 288c 289 if(me.eq.0) then 290c 291 open(unit=lfnrst,file=filrst(1:index(filrst,' ')-1), 292 + status='old',form='formatted',err=9999) 293 rewind(lfnrst) 294c 295 do 2 i=1,3 296 read(lfnrst,1001) cdum 297 1001 format(a1) 298 2 continue 299 read(lfnrst,1006) nhist 300 1006 format(32x,i5) 301 do 6 i=1,nhist 302 read(lfnrst,1007) hist(i) 303 1007 format(a) 304 6 continue 305 read(lfnrst,1002) npbtyp,nbxtyp,rsgmr 306 1002 format(2i5,f12.6) 307 if(rsgm.lt.0.0d0) rsgm=rsgmr 308 read(lfnrst,1004) ((vlat(i,j),j=1,3),i=1,3) 309 1004 format(3f12.6) 310 box(1)=vlat(1,1) 311 box(2)=vlat(2,2) 312 box(3)=vlat(3,3) 313 read(lfnrst,1003) jdum 314 1003 format(40x,i5) 315 read(lfnrst,1001) cdum 316 if(jdum.ne.0) then 317 read(lfnrst,1001) cdum 318 endif 319 read(lfnrst,1005) nwm,nwa,nsm,nsa,nwmc,nsf,nseq 320 1005 format(7i10) 321 close(unit=lfnrst,status='keep') 322c 323 bsize=max(rshort+half*rsgm,half*(rlong+half*rsgm),rbox) 324c 325 if(util_print('distribution',print_default)) then 326 if(me.eq.0) write(lfnout,2001) 327 2001 format(/,' Distribution information',/) 328 if(me.eq.0) write(lfnout,2002) rshort,rsgm,rlong,rbox,box,bsize 329 2002 format(' Short range cutoff',t35,f12.6,/, 330 + ' Segment size',t35,f12.6,/,' Long range cutoff',t35,f12.6,/, 331 + ' Box size rbox',t35,f12.6,//,' Box dimension',t35,3f12.6,//, 332 + ' Initial cell size',t35,f12.6,/) 333 endif 334c 335 if(nbx*nby*nbz.lt.np) then 336 nbx=int(box(1)/bsize) 337 nby=int(box(2)/bsize) 338 nbz=int(box(3)/bsize) 339 endif 340c 341 nbx=max(1,nbx,npx) 342 nby=max(1,nby,npy) 343 nbz=max(1,nbz,npz) 344c 345 if(util_print('distribution',print_default)) then 346 if(me.eq.0) then 347 write(lfnout,2003) nbx,nby,nbz 348 2003 format(' Initial cell distribution',t35,3i5) 349 endif 350 endif 351c 352 if(nbxin.eq.0) then 353 nred(1)=0 354 nred(2)=0 355 nred(3)=0 356 if(nbx.gt.npx.and.mod(nbx,npx).gt.0) then 357 nbx=(nbx/npx)*npx 358 nred(1)=nbx 359 endif 360 endif 361 if(nbyin.eq.0) then 362 if(nby.gt.npy.and.mod(nby,npy).gt.0) then 363 nby=(nby/npy)*npy 364 nred(2)=nby 365 endif 366 endif 367 if(nbzin.eq.0) then 368 if(nbz.gt.npz.and.mod(nbz,npz).gt.0) then 369 nbz=(nbz/npz)*npz 370 nred(3)=nbz 371 endif 372 endif 373c 374 if(util_print('distribution',print_default)) then 375 if(me.eq.0) then 376 write(lfnout,2004) nbx,nby,nbz 377 2004 format(' Final cell distribution',t35,3i5,/) 378 endif 379 endif 380c 381 bxmin=bsize/dble(int((dble(nbx)*bsize)/box(1))+1) 382 bymin=bsize/dble(int((dble(nby)*bsize)/box(2))+1) 383 bzmin=bsize/dble(int((dble(nbz)*bsize)/box(3))+1) 384c 385 if(util_print('distribution',print_default)) then 386 if(me.eq.0) then 387 write(lfnout,2005) bxmin,bymin,bzmin 388 2005 format(' Minimum cell size',t35,3f12.6) 389 endif 390 endif 391c 392 endif 393c 394 if(np.gt.1) then 395 itemp(1)=nwm 396 itemp(2)=nwa 397 itemp(3)=nsm 398 itemp(4)=nsa 399 itemp(5)=nbx 400 itemp(6)=nby 401 itemp(7)=nbz 402 itemp(8)=npbtyp 403 itemp(9)=nbxtyp 404 itemp(10)=nsf 405 itemp(11)=nseq 406 rtemp(1)=bsize 407 rtemp(2)=1.001d0*bxmin 408 rtemp(3)=1.001d0*bymin 409 rtemp(4)=1.001d0*bzmin 410 rtemp(5)=rsgm 411 call ga_brdcst(msp_03,itemp,ma_sizeof(mt_int,11,mt_byte),0) 412 call ga_brdcst(msp_04,rtemp,ma_sizeof(mt_dbl,5,mt_byte),0) 413 call ga_brdcst(msp_05,box,ma_sizeof(mt_dbl,3,mt_byte),0) 414 call ga_brdcst(msp_06,vlat,ma_sizeof(mt_dbl,9,mt_byte),0) 415 nwm=itemp(1) 416 nwa=itemp(2) 417 nsm=itemp(3) 418 nsa=itemp(4) 419 nbx=itemp(5) 420 nby=itemp(6) 421 nbz=itemp(7) 422 npbtyp=itemp(8) 423 nbxtyp=itemp(9) 424 nsf=itemp(10) 425 nseq=itemp(11) 426 bsize=rtemp(1) 427 bxmin=rtemp(2) 428 bymin=rtemp(3) 429 bzmin=rtemp(4) 430 rsgm=rtemp(5) 431 endif 432c 433 rbbl=rlong+half*rsgm 434c 435 boxh(1)=half*box(1) 436 boxh(2)=half*box(2) 437 boxh(3)=half*box(3) 438c 439 maxbox=max(nbx,nby,nbz) 440 nbtot=nbx*nby*nbz 441 lpbc0=nbx.eq.1.or.nby.eq.1.or.nbz.eq.1.or. 442 + npx.eq.1.or.npy.eq.1.or.npz.eq.1.or.lpbc9 443c 444 mbox=30 445 do 3 ibx=1,nbx 446 do 4 iby=1,nby 447 do 5 ibz=1,nbz 448 if(me.eq.((((ibz-1)*npz)/nbz)*npy+(((iby-1)*npy)/nby))*npx 449 + +((ibx-1)*npx)/nbx) mbox=mbox+1 450 5 continue 451 4 continue 452 3 continue 453 mbxloc=mbox-30 454c 455 if(np.gt.1) call ga_igop(msp_07,mbox,1,'max') 456c 457 mwa=max(1,nwa) 458 msag=max(1,msa,(mbox-30+madbox)*((nwm*nwa+nsa)/nbtot+1)+1) 459 msa=max(1,msa,(mbox-30+madbox)*((nwm*nwa+nsa)/nbtot+1)+1) 460 mwmg=max(1,msag/mwa+1) 461 mwm=max(1,mwm,msa/mwa+1) 462c 463 msa=min(msa,2*nsa+1) 464 mwm=min(mwm,2*nwm+1) 465 msag=min(msag,2*nsa+1) 466 mwmg=min(mwmg,2*nwm+1) 467c 468 if(lnode0) then 469 msa=nsa+1 470 mwm=nwm+1 471 msag=nsa+1 472 mwmg=nwm+1 473 endif 474c 475 if(util_print('distribution',print_default)) then 476 if(me.eq.0) then 477 write(lfnout,2006) mbox-30,madbox,nbtot 478 2006 format(/,' ARRAY DIMENSION INFORMATION',//, 479 + ' Number cells per processor: ',i7,/, 480 + ' Number of buffer cells: ',i7,/, 481 + ' Total number of cells: ',i7) 482 endif 483 endif 484c 485 return 486 9999 call md_abort('Failed to open restart file',0) 487 return 488 end 489 subroutine sp_alloc 490c 491 implicit none 492c 493#include "sp_common.fh" 494#include "mafdecls.fh" 495#include "global.fh" 496c 497 integer isize 498c 499 if(.not.ma_push_get(mt_int,np,'bindex',l_bindex,i_bindex)) 500 + call md_abort('Failed to allocate bindex',0) 501 if(.not.ma_push_get(mt_int,54*np,'buren',l_buren,i_buren)) 502 + call md_abort('Failed to allocate buren',0) 503 if(.not.ma_push_get(mt_int,3*maxbox,'owner',l_iown,i_iown)) 504 + call md_abort('Failed to allocate owner',0) 505 if(.not.ma_push_get(mt_dbl,3*maxbox,'bxsiz',l_boxs,i_boxs)) 506 + call md_abort('Failed to allocate bxsiz',0) 507 if(.not.ma_push_get(mt_int,6*maxbox,'ibxrg',l_boxr,i_boxr)) 508 + call md_abort('Failed to allocate ibxrg',0) 509 if(.not.ma_push_get(mt_int,6*maxbox,'rng',l_rng,i_rng)) 510 + call md_abort('Failed to allocate rng',me) 511c 512 if(.not.ma_push_get(mt_int,mip2*mbox,'ipl',l_ipl,i_ipl)) 513 + call md_abort('Failed to allocate ipl',0) 514 if(.not.ma_push_get(mt_int,mip2*mbox,'jpl',l_jpl,i_jpl)) 515 + call md_abort('Failed to allocate jpl',0) 516c 517 mseq=nseq 518c 519 if(.not.ma_push_get(mt_int,mseq,'lseq',l_lseq,i_lseq)) 520 + call md_abort('Failed to allocate lseq',0) 521 if(.not.ma_push_get(mt_int,mseq,'sndx',l_sndx,i_sndx)) 522 + call md_abort('Failed to allocate sndx',0) 523c 524 if(.not.ga_create(mt_int,np*mbox,mip2,'ip',mbox,mip2,ga_ip)) 525 + call md_abort('Failed to create global array ip',0) 526c 527 return 528 end 529 subroutine sp_alloc2 530c 531 implicit none 532c 533#include "sp_common.fh" 534#include "mafdecls.fh" 535#include "global.fh" 536#include "util.fh" 537c 538 integer isize 539c 540 if(nbget.ne.0) then 541 msa=max(1,msa,(mbox-30+nrempr+madbox)*((nwm*nwa+nsa)/nbtot+1)+1) 542 mwm=max(1,mwm,msa/mwa+1) 543 msa=min(msa,2*nsa+1) 544 mwm=min(mwm,2*nwm+1) 545 endif 546c 547 if(util_print('distribution',print_default)) then 548 if(me.eq.0) then 549 if(nbget.ne.0) then 550 write(lfnout,2005) nrempr 551 2005 format(' Number of remote cell pairs: ',i7) 552 if(nbget.gt.0) then 553 write(lfnout,2006) nbget 554 2006 format(' Number of prefetch cells: ',i7) 555 endif 556 endif 557 write(lfnout,2007) mwm,mwmg 558 2007 format(' Dimension solvent local: ',i7,', global:',i7) 559 write(lfnout,2008) msa,msag 560 2008 format(' Dimension solute local: ',i7,', global:',i7) 561 endif 562 endif 563c 564 if(.not.ga_create(mt_int,np*mwmg,miw2,'iw',mwmg,miw2,ga_iw)) 565 + call md_abort('Failed to create global array iw',0) 566 isize=6+12*mwa 567 if(lpola) isize=6+18*mwa 568 if(lpola.and.lfree) isize=6+30*mwa 569 if(.not.ga_create(mt_dbl,np*mwmg,isize,'w',mwmg,isize,ga_w)) 570 + call md_abort('Failed to create global array w',0) 571 if(.not.ga_create(mt_int,np*msag,mis2,'is',msag,mis2,ga_is)) 572 + call md_abort('Failed to create global array is',0) 573 isize=39 574 if(lpola) isize=45 575 if(lpola.and.lfree) isize=57 576 if(.not.ga_create(mt_dbl,np*msag,isize,'s',msag,isize,ga_s)) 577 + call md_abort('Failed to create global array s',0) 578c 579 if(.not.ga_create(mt_int,np*mwmg,1,'iwz',mwmg,1,ga_iwz)) 580 + call md_abort('Failed to create global array iwz',0) 581 if(.not.ga_create(mt_int,np*msag,1,'isz',msag,1,ga_isz)) 582 + call md_abort('Failed to create global array isz',0) 583c 584 return 585 end 586 subroutine sp_free() 587c 588 implicit none 589c 590#include "sp_common.fh" 591#include "mafdecls.fh" 592#include "global.fh" 593c 594 if(.not.ga_destroy(ga_isz)) 595 + call md_abort('Error ga_destroy isz',ga_isz) 596 if(.not.ga_destroy(ga_iwz)) 597 + call md_abort('Error ga_destroy iwz',ga_iwz) 598c 599 if(.not.ga_destroy(ga_s)) 600 + call md_abort('Error ga_destroy s',ga_s) 601 if(.not.ga_destroy(ga_is)) 602 + call md_abort('Error ga_destroy is',ga_is) 603 if(.not.ga_destroy(ga_w)) 604 + call md_abort('Error ga_destroy w',ga_w) 605 if(.not.ga_destroy(ga_iw)) 606 + call md_abort('Error ga_destroy iw',ga_iw) 607 if(.not.ga_destroy(ga_ip)) 608 + call md_abort('Error ga_destroy iw',ga_ip) 609c 610 if(.not.ma_pop_stack(l_sndx)) 611 + call md_abort('Failed to deallocate sndx',0) 612 if(.not.ma_pop_stack(l_lseq)) 613 + call md_abort('Failed to deallocate lseq',0) 614 if(.not.ma_pop_stack(l_jpl)) 615 + call md_abort('Failed to deallocate jpl',0) 616 if(.not.ma_pop_stack(l_ipl)) 617 + call md_abort('Failed to deallocate ipl',0) 618 if(.not.ma_pop_stack(l_rng)) 619 + call md_abort('Failed to deallocate rng',0) 620 if(.not.ma_pop_stack(l_boxr)) 621 + call md_abort('Failed to deallocate boxr',0) 622 if(.not.ma_pop_stack(l_boxs)) 623 + call md_abort('Failed to deallocate boxs',0) 624 if(.not.ma_pop_stack(l_iown)) 625 + call md_abort('Failed to deallocate iown',0) 626 if(.not.ma_pop_stack(l_buren)) 627 + call md_abort('Failed to deallocate buren',0) 628 if(.not.ma_pop_stack(l_bindex)) 629 + call md_abort('Failed to deallocate bindex',0) 630c 631 return 632 end 633 subroutine sp_rdrst(lfnrst,filrst,lfntop,filtop, 634 + temp,tempw,temps,ipl,xw,vw,fw,xwcr,iwl,xs,vs,fs,xscr,isl, 635 + bxw,bvw,bfw,brw,ibw,nw,bxs,bvs,bfs,ibs,ns,ndx, 636 + ibownr,boxsiz,lseq,isndx) 637c 638 implicit none 639c 640#include "sp_common.fh" 641#include "mafdecls.fh" 642#include "msgids.fh" 643c 644 integer sp_btop 645 external sp_btop 646c 647 integer lfnrst,lfntop 648 character*255 filrst,filtop 649 integer nw,ns 650 integer ibownr(maxbox,3) 651 integer ipl(mbox,mip2),iwl(mwm,miw2),isl(msa,mis2) 652 real*8 xw(mwm,3,mwa),xs(msa,3),xwcr(mwm,3) 653 real*8 vw(mwm,3,mwa),vs(msa,3),xscr(msm,3) 654 real*8 fw(mwm,3,mwa),fs(msa,3) 655 real*8 bxw(nw,3,mwa),bxs(ns,3),brw(nw,3) 656 real*8 bvw(nw,3,mwa),bvs(ns,3) 657 real*8 bfw(nw,3,mwa),bfs(ns,3) 658 integer ibw(nw,2),ibs(ns,mis2) 659 real*8 boxsiz(maxbox,3) 660 integer ndx(nw),lseq(mseq),isndx(mseq) 661 real*8 temp,tempw,temps 662c 663 character*1 cdum 664 real*8 rdum,cgx,cgy,cgz 665 integer i,j,k,idum,jdum,kdum,number,ncyc,numw 666 integer icyc,ibx,iby,ibz,ipx,ipy,ipz,node,new,nold 667 integer ilw,ihw,jlw,jhw,ils,ihs,jls,jhs 668 integer ili,ihi,jli,jhi,ilp,ihp,jlp,jhp 669c 670 integer nat,nqt,naw,nbw,nhw,ndw,now,ntw,nnw,nsmr 671 real*8 boxi(3) 672 logical lforces 673 character*80 card 674c 675 integer l,m,ib(3),nbox,nb(3),joff,npars,icount,lasts 676 real*8 xtmin,xtmax,xtx,xt(3) 677c 678 if(.not.ma_verify_allocator_stuff()) then 679 call md_abort('ERROR IN MEMORY USE ENTERING SP_RDRST',0) 680 endif 681c 682 lpbc9=.true. 683 nsmr=0 684 icount=0 685 lasts=0 686 ndums=0 687c 688 if(me.eq.0) then 689c 690 open(unit=lfnrst,file=filrst(1:index(filrst,' ')-1), 691 + status='old',form='formatted',err=9999) 692 rewind(lfnrst) 693c 694 open(unit=lfntop,file=filtop(1:index(filtop,' ')-1), 695 + status='old',form='formatted',err=9899) 696 rewind(lfntop) 697c 698 write(lfnout,6000) 699 6000 format(/,' TOPOLOGY FILE INFORMATION',/) 700c 701 read(lfntop,2001,end=9897,err=9898) card 702 2001 format(a80) 703 write(lfnout,6001) card 704 6001 format(' Title',t15,a) 705 read(lfntop,2001,end=9897,err=9898) card 706 write(lfnout,6002) card 707 6002 format(t15,a) 708 read(lfntop,2001,end=9897,err=9898) card 709 write(lfnout,6002) card 710c 711 read(lfntop,2001,end=9897,err=9898) card 712 write(lfnout,6003) card(1:12),card(13:32),card(33:42) 713 6003 format(' Version',t15,a,/,' Date',t15,a,/,' Force field',t15,a) 714c 715 read(lfntop,2002,end=9897,err=9898) npars 716 read(lfntop,2002,end=9897,err=9898) nat 717 read(lfntop,2002,end=9897,err=9898) nqt 718 read(lfntop,2002,end=9897,err=9898) nseq 719 2002 format(i5) 720 read(lfntop,2001,end=9897,err=9898) 721 do 1103 k=1,npars 722 do 102 i=1,nat 723 read(lfntop,2001,end=9897,err=9898) 724 do 103 j=i,nat 725 read(lfntop,2001,end=9897,err=9898) 726 103 continue 727 102 continue 728 1103 continue 729 do 104 i=1,nqt*npars 730 read(lfntop,2001,end=9897,err=9898) 731 104 continue 732 do 4104 i=1,nseq 733 read(lfntop,2001,end=9897,err=9898) 734 4104 continue 735 read(lfntop,2003,end=9897,err=9898) naw,nbw,nhw,ndw,now,ntw,nnw 736 2003 format(5i7,2i10) 737 read(lfntop,2001,end=9897,err=9898) 738 do 105 i=1,naw 739 read(lfntop,2001,end=9897,err=9898) 740 105 continue 741 do 106 i=1,nbw*(npars+1) 742 read(lfntop,2001,end=9897,err=9898) 743 106 continue 744 do 107 i=1,nhw*(npars+1) 745 read(lfntop,2001,end=9897,err=9898) 746 107 continue 747 do 108 i=1,ndw*(npars+1) 748 read(lfntop,2001,end=9897,err=9898) 749 108 continue 750 do 109 i=1,now*(npars+1) 751 read(lfntop,2001,end=9897,err=9898) 752 109 continue 753 if(ntw.gt.0) then 754 read(lfntop,2004,end=9897,err=9898) 755 read(lfntop,2004,end=9897,err=9898) 756 2004 format(11i7) 757 endif 758 if(nnw.gt.0) then 759 read(lfntop,2005,end=9997,err=9998) 760 read(lfntop,2005,end=9997,err=9998) 761 2005 format(11i7) 762 endif 763 read(lfntop,2001,end=9897,err=9898) 764 do 204 i=1,npars 765 read(lfntop,2001,end=9897,err=9898) 766 204 continue 767c 768 write(lfnout,6100) 769 6100 format(/,' RESTART FILE INFORMATION',/) 770c 771 read(lfnrst,1001,end=9997,err=9998) card 772 1001 format(a80) 773 write(lfnout,6101) card 774 6101 format(' Title',t15,a) 775 read(lfnrst,1001,end=9997,err=9998) card 776 write(lfnout,6102) card 777 6102 format(t15,a) 778 read(lfnrst,1001,end=9997,err=9998) card 779 write(lfnout,6002) card 780c 781 read(lfnrst,1016) card(1:32),nhist,lforces 782 1016 format(a32,i5,4x,l1) 783 write(lfnout,6103) card(1:12),card(13:32) 784 6103 format(' Version',t15,a,/,' Date',t15,a) 785c 786 if(nhist.gt.0) write(lfnout,6104) 787 6104 format(/,' History',/) 788c 789 if(nhist.gt.0) then 790 do 21 i=1,nhist 791 read(lfnrst,1017) hist(i) 792 1017 format(a) 793 write(lfnout,6105) hist(i) 794 6105 format(1x,a) 795 21 continue 796 endif 797 if(nhist.lt.mxhist) then 798 nhist=nhist+1 799 else 800 do 22 i=1,nhist-1 801 hist(i)=hist(i+1) 802 22 continue 803 endif 804 do 23 i=1,80 805 hist(nhist)(i:i)=' ' 806 23 continue 807 read(lfnrst,1002,end=9997,err=9998) npbtyp,nbxtyp 808 1002 format(i5,i5) 809 read(lfnrst,1003,end=9997,err=9998) ((vlat(i,j),j=1,3),i=1,3) 810 1003 format(3f12.6) 811 read(lfnrst,1004,end=9997,err=9998) jdum 812 1004 format(40x,i5) 813 read(lfnrst,1005,end=9997,err=9998) temp,tempw,temps 814 1005 format(3f12.6) 815 do 2 i=1,3 816 box(i)=vlat(i,i) 817 boxh(i)=half*box(i) 818 boxi(i)=one/box(i) 819 do 3 j=1,3 820 vlati(i,j)=vlat(i,j) 821 3 continue 822 2 continue 823 call matinv(vlati,3,3) 824 if(jdum.ne.0) then 825 read(lfnrst,1001,end=9997,err=9998) cdum 826 endif 827 read(lfnrst,1006,end=9997,err=9998) idum 828 1006 format(70x,i5) 829 if(idum.gt.0) then 830 read(lfnrst,1007,end=9997,err=9998) idum,jdum,kdum 831 1007 format(3i5) 832 read(lfnrst,1008,end=9997,err=9998) (rdum,i=1,idum) 833 read(lfnrst,1008,end=9997,err=9998) (rdum,i=1,jdum) 834 read(lfnrst,1008,end=9997,err=9998) (rdum,i=1,kdum) 835 1008 format(4e20.12) 836 endif 837c 838 if(nwm.gt.0) then 839 number=0 840 ncyc=nwm/nw+1 841 numw=nw 842 do 4 icyc=1,ncyc 843 if(nwm-number.lt.numw) numw=nwm-number 844 do 44 i=1,numw 845 read(lfnrst,1009,end=9997,err=9998) 846 + ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa) 847 1009 format(2x,6f13.8) 848 if(lforces) read(lfnrst,1109,end=9997,err=9998) 849 + ((bfw(i,j,k),j=1,3),k=1,nwa) 850 1109 format(2x,6e13.6) 851 read(lfnrst,1010,end=9997,err=9998) ibw(i,2),(brw(i,k),k=1,3) 852 1010 format(i1,1x,3f13.8) 853 44 continue 854 do 5 i=1,numw 855 cgx=zero 856 cgy=zero 857 cgz=zero 858 do 6 k=1,nwa 859 cgx=cgx+bxw(i,1,k) 860 cgy=cgy+bxw(i,2,k) 861 cgz=cgz+bxw(i,3,k) 862 if(.not.lforces) then 863 bfw(i,1,k)=zero 864 bfw(i,2,k)=zero 865 bfw(i,3,k)=zero 866 endif 867 6 continue 868 ibx=0 869 iby=0 870 ibz=0 871 if(nbxtyp.ne.1) then 872 xt(1)=cgx 873 xt(2)=cgy 874 xt(3)=cgz 875 else 876 xt(1)=box(1)*(vlati(1,1)*cgx+vlati(1,2)*cgy+vlati(1,3)*cgz) 877 xt(2)=box(2)*(vlati(2,1)*cgx+vlati(2,2)*cgy+vlati(2,3)*cgz) 878 xt(3)=box(3)*(vlati(3,1)*cgx+vlati(3,2)*cgy+vlati(3,3)*cgz) 879 endif 880 do 7 j=1,nbx-1 881 if(xt(1)/nwa+boxh(1).gt.boxsiz(j,1)) ibx=j 882 7 continue 883 do 8 j=1,nby-1 884 if(xt(2)/nwa+boxh(2).gt.boxsiz(j,2)) iby=j 885 8 continue 886 do 9 j=1,nbz-1 887 if(xt(3)/nwa+boxh(3).gt.boxsiz(j,3)) ibz=j 888 9 continue 889c 890 if(npbtyp.gt.0) then 891 if(ibx.ge.nbx) ibx=ibx-nbx 892 if(iby.ge.nby) iby=iby-nby 893 if(ibx.lt.0) ibx=ibx+nbx 894 if(iby.lt.0) iby=iby+nby 895 if(npbtyp.eq.1) then 896 if(ibz.ge.nbz) ibz=ibz-nbz 897 if(ibz.lt.0) ibz=ibz+nbz 898 else 899 if(ibz.ge.nbz) ibz=nbz-1 900 if(ibz.lt.0) ibz=0 901 endif 902 else 903 if(ibx.ge.nbx) ibx=nbx-1 904 if(iby.ge.nby) iby=nby-1 905 if(ibz.ge.nbz) ibz=nbz-1 906 if(ibx.lt.0) ibx=0 907 if(iby.lt.0) iby=0 908 if(ibz.lt.0) ibz=0 909 endif 910 ipx=ibownr(ibx+1,1) 911 ipy=ibownr(iby+1,2) 912 ipz=ibownr(ibz+1,3) 913c 914 ndx(i)=(ipz*npy+ipy)*npx+ipx 915 ibw(i,1)=(ibz*nby+iby)*nbx+ibx 916 if(lnode0) then 917 ndx(i)=0 918 ibw(i,1)=0 919 endif 920 5 continue 921 do 10 node=0,np-1 922 new=0 923 do 11 i=1,numw 924 if(ndx(i).eq.node) then 925 new=new+1 926 iwl(new,lwgmn)=number+i 927 iwl(new,lwnod)=node 928 do 124 j=1,3 929 do 12 k=1,nwa 930 xw(new,j,k)=bxw(i,j,k) 931 vw(new,j,k)=bvw(i,j,k) 932 fw(new,j,k)=bfw(i,j,k) 933 12 continue 934 if(nserie.eq.0) then 935 xwcr(new,j)=zero 936 else 937 xwcr(new,j)=brw(i,j) 938 endif 939 124 continue 940 iwl(new,lwbox)=ibw(i,1) 941 iwl(new,lwdyn)=5*ibw(i,2) 942 endif 943 11 continue 944c 945 if(new.gt.0) then 946 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 947 call ga_get(ga_ip,ilp,ilp,jlp,jhp,ipl,mbox) 948 nold=ipl(1,2) 949 if(nold+new.gt.mwm) call md_abort('Dimension mwm too small',0) 950 ipl(1,2)=ipl(1,2)+new 951 call ga_put(ga_ip,ilp,ilp,jlp,jhp,ipl,mbox) 952 call ga_distribution(ga_iw,node,ili,ihi,jli,jhi) 953 ili=ili+nold 954 ihi=ili+new-1 955 call ga_put(ga_iw,ili,ihi,jli,jhi,iwl,mwm) 956 call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw) 957 ilw=ilw+nold 958 ihw=ilw+new-1 959 call ga_put(ga_w,ilw,ihw,jlw,jlw+3*mwa-1,xw(1,1,1),mwm) 960 call ga_put(ga_w,ilw,ihw,jlw+3*mwa,jlw+6*mwa-1,vw(1,1,1),mwm) 961 call ga_put(ga_w,ilw,ihw,jlw+6*mwa,jlw+6*mwa+2,xwcr(1,1),mwm) 962 call ga_put(ga_w,ilw,ihw,jlw+6*mwa+3,jlw+9*mwa+2,fw(1,1,1),mwm) 963 endif 964c 965 10 continue 966 number=number+numw 967 4 continue 968 endif 969c 970 if(nsa.gt.0) then 971 nb(1)=nbx 972 nb(2)=nby 973 nb(3)=nbz 974 nbox=0 975 k=1 976 joff=0 977 do 13 i=1,nsa 978 if(k.gt.ns) call md_abort('Increase memory for sp_rdrst buffer',0) 979 read(lfnrst,1011,end=9997,err=9998) 980 + ibs(k,11),(bxs(k,j),j=1,3),(bvs(k,j),j=1,3),ibs(k,12) 981 1011 format(i1,1x,6f13.8,i5) 982 if(ibs(k,12).lt.0) ndums=ndums+1 983 if(lforces) then 984 read(lfnrst,1111,end=9997,err=9998) (bfs(k,j),j=1,3) 985 1111 format(2x,6e13.6) 986 else 987 bfs(k,1)=zero 988 bfs(k,2)=zero 989 bfs(k,3)=zero 990 endif 991 read(lfntop,2009,end=9897,err=9898) (ibs(k,j),j=1,10) 992 2009 format(16x,i3,4i7,5i5) 993c 2009 format(16x,10i5) 994 if(ibs(k,3).ne.lasts) then 995 icount=icount+1 996 lasts=ibs(k,3) 997 isndx(icount)=lasts 998 endif 999 ibs(k,3)=icount 1000 if(nsmr.lt.ibs(k,2)) nsmr=ibs(k,2) 1001 if(nserie.eq.0) then 1002 if(nsf.lt.ibs(k,1)) nsf=ibs(k,1) 1003 else 1004 if(nsf.lt.ibs(k,1)) 1005 + call md_abort('Error in number of solute fractions',nsf) 1006 endif 1007c 1008c if segment of this atom is different from the segment of previous atom 1009c then distribute all previous atoms in the list 1010c 1011 if(k.gt.1) then 1012 if(ibs(k,3).ne.ibs(k-1,3)) then 1013 new=k-1 1014 goto 14 1015 endif 1016 endif 1017c 1018c if this is the last atom distribute 1019c 1020 if(i.eq.nsa) then 1021 new=k 1022 goto 14 1023 endif 1024c 1025c read next atom 1026c 1027 k=k+1 1028 goto 13 1029c 1030c distribute atoms 1 through new 1031c 1032 14 continue 1033c 1034c determine the center of geometry 1035c 1036 do 15 l=1,3 1037 xtmax=bxs(1,l) 1038 xtmin=bxs(1,l) 1039 do 16 j=1,new 1040 xtmax=max(xtmax,bxs(j,l)) 1041 xtmin=min(xtmin,bxs(j,l)) 1042 16 continue 1043 xtx=0.5d0*(xtmax+xtmin) 1044 if(npbtyp.ne.0) then 1045 if(abs(xtx).gt.boxh(l)) then 1046 xtx=xtx-nint(xtx*boxi(l))*box(l) 1047 endif 1048 endif 1049 ib(l)=0 1050 do 17 m=1,nb(l)-1 1051 if(xtx+boxh(l).gt.boxsiz(m,l)) ib(l)=m 1052 17 continue 1053 15 continue 1054c 1055c periodic boundaries 1056c 1057 if(npbtyp.gt.0) then 1058 m=2 1059 if(npbtyp.eq.1) m=3 1060 do 18 l=1,m 1061 if(ib(l).ge.nb(l)) ib(l)=ib(l)-nb(l) 1062 if(ib(l).lt.0) ib(l)=ib(l)+nb(l) 1063 18 continue 1064 if(npbtyp.gt.1) then 1065 if(ib(3).ge.nb(3)) ib(3)=nb(3)-1 1066 if(ib(3).lt.0) ib(3)=0 1067 endif 1068 else 1069 do 19 l=1,3 1070 if(ib(l).ge.nb(l)) ib(l)=nb(l)-1 1071 if(ib(l).lt.0) ib(l)=0 1072 19 continue 1073 endif 1074c 1075c determine owning node 1076c 1077 if(.not.lnode0) nbox=(ib(3)*nb(2)+ib(2))*nb(1)+ib(1) 1078 node=sp_btop(nbox,ibownr) 1079c 1080 do 120 j=1,new 1081 isl(j,lsgan)=joff+j 1082 isl(j,lsfrc)=ibs(j,1) 1083 isl(j,lsmol)=ibs(j,2) 1084 isl(j,lssgm)=ibs(j,3) 1085 isl(j,lsgrp)=ibs(j,4) 1086 isl(j,lspgr)=ibs(j,5) 1087 isl(j,lsatt)=ibs(j,6) 1088 isl(j,lsct1)=ibs(j,7) 1089 isl(j,lsct2)=ibs(j,8) 1090 isl(j,lsct3)=ibs(j,9) 1091 isl(j,lssss)=ibs(j,10) 1092 isl(j,lsdyn)=5*ibs(j,11) 1093 isl(j,lsbox)=nbox 1094 isl(j,lsnod)=node 1095 isl(j,lshop)=0 1096 if(ibs(j,12).gt.0) isl(j,lshop)=ibs(j,12)*2 1097 if(ibs(j,12).lt.0) isl(j,lshop)=(-ibs(j,12))*2+1 1098 do 121 l=1,3 1099 xs(j,l)=bxs(j,l) 1100 vs(j,l)=bvs(j,l) 1101 fs(j,l)=bfs(j,l) 1102 121 continue 1103 120 continue 1104 joff=joff+new 1105c 1106c communicate data to node 1107c 1108 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 1109 call ga_get(ga_ip,ilp,ilp+1,jlp,jhp,ipl,mbox) 1110 nold=ipl(2,2) 1111 if(nold+new.gt.msa) 1112 + call md_abort('Dimension msa too small (1)',nold+new) 1113 ipl(2,2)=ipl(2,2)+new 1114 call ga_put(ga_ip,ilp,ilp+1,jlp,jhp,ipl,mbox) 1115 call ga_distribution(ga_is,node,ili,ihi,jli,jhi) 1116 ili=ili+nold 1117 ihi=ili+new-1 1118 call ga_put(ga_is,ili,ihi,jli,jhi,isl,msa) 1119 call ga_distribution(ga_s,node,ils,ihs,jls,jhs) 1120 ils=ils+nold 1121 ihs=ils+new-1 1122 call ga_put(ga_s,ils,ihs,jls,jls+2,xs(1,1),msa) 1123 call ga_put(ga_s,ils,ihs,jls+3,jls+5,vs(1,1),msa) 1124 call ga_put(ga_s,ils,ihs,jls+9,jls+11,fs(1,1),msa) 1125c 1126c make first atom of next segment first in the list 1127c 1128 if(k.gt.new) then 1129 do 122 j=1,12 1130 ibs(1,j)=ibs(k,j) 1131 122 continue 1132 do 123 j=1,3 1133 bxs(1,j)=bxs(k,j) 1134 bvs(1,j)=bvs(k,j) 1135 bfs(1,j)=bfs(k,j) 1136 123 continue 1137 k=2 1138 if(i.eq.nsa) then 1139 new=1 1140 k=1 1141 goto 14 1142 endif 1143 endif 1144c 1145 13 continue 1146 endif 1147c 1148 if(nsm.gt.0) then 1149 do 31 i=1,nsm 1150 xscr(i,1)=zero 1151 xscr(i,2)=zero 1152 xscr(i,3)=zero 1153 31 continue 1154 do 32 i=1,nsm 1155 read(lfnrst,1012,end=99,err=99) (xscr(i,j),j=1,3) 1156 1012 format(2x,3f13.8) 1157 if(nserie.eq.0) then 1158 xscr(i,1)=zero 1159 xscr(i,2)=zero 1160 xscr(i,3)=zero 1161 endif 1162 32 continue 1163 endif 1164c 1165 if(.not.ma_verify_allocator_stuff()) then 1166 call md_abort('ERROR IN MEMORY IN SP_RDRST',1) 1167 endif 1168c 1169 if(nseq.gt.0) then 1170 read(lfnrst,1013) (lseq(i),i=1,nseq) 1171 1013 format(20i3) 1172 endif 1173c 1174 if(.not.ma_verify_allocator_stuff()) then 1175 call md_abort('ERROR IN MEMORY IN SP_RDRST',2) 1176 endif 1177c 1178 99 continue 1179c 1180 close(unit=lfnrst,status='keep') 1181 close(unit=lfntop,status='keep') 1182c 1183 endif 1184c 1185 if(nseq.gt.0) then 1186 call ga_brdcst(msp_22,lseq,nseq*ma_sizeof(mt_int,1,mt_byte),0) 1187 endif 1188 call ga_brdcst(msp_28,ndums,ma_sizeof(mt_int,1,mt_byte),0) 1189c 1190 if(.not.ma_verify_allocator_stuff()) then 1191 call md_abort('ERROR IN MEMORY USE EXITING SP_RDRST',0) 1192 endif 1193c 1194 return 1195 9897 call md_abort('EOF encountered on topology file',0) 1196 9898 call md_abort('Error reading topology file',1) 1197 9899 call md_abort('Error opening topology file',2) 1198 9997 call md_abort('EOF encountered on restart file',0) 1199 9998 call md_abort('Error reading restart file',0) 1200 9999 call md_abort('Error opening restart file',0) 1201 return 1202 end 1203 integer function sp_btop(ibox,ibownr) 1204c 1205 implicit none 1206c 1207#include "sp_common.fh" 1208c 1209 integer ibox,ibownr(maxbox,3) 1210 integer iboxx,iboxy,iboxz 1211c 1212 iboxx=mod(ibox,nbx) 1213 iboxy=mod((ibox-iboxx)/nbx,nby) 1214 iboxz=((ibox-iboxx)/nbx-iboxy)/nby 1215 sp_btop=(ibownr(iboxz+1,3)*npy+ibownr(iboxy+1,2))*npx+ 1216 + ibownr(iboxx+1,1) 1217c 1218 return 1219 end 1220 subroutine sp_decomp(ibownr,boxsiz,iburen,ibindx) 1221c 1222 implicit none 1223c 1224#include "sp_common.fh" 1225#include "msgids.fh" 1226#include "global.fh" 1227c 1228 integer ibownr(maxbox,3) 1229 real*8 boxsiz(maxbox,3) 1230 integer iburen(np,27,2) 1231 integer ibindx(np) 1232c 1233 integer ibx,iby,ibz,i,j,ix,jx,iy,jy,iz,jz,jnode,nrnod,nbiown 1234c 1235c check dimensions of ibownr 1236c 1237 if(nbx.gt.maxbox.or.nby.gt.maxbox.or.nbz.gt.maxbox) 1238 + call md_abort('Dimension maxbox too small',0) 1239c 1240c determine the node dimension for each sub box 1241c 1242 do 1 ibx=1,nbx 1243 ibownr(ibx,1)=((ibx-1)*npx)/nbx 1244 boxsiz(ibx,1)=(box(1)*dble(ibx))/dble(nbx) 1245 1 continue 1246 do 2 iby=1,nby 1247 ibownr(iby,2)=((iby-1)*npy)/nby 1248 boxsiz(iby,2)=(box(2)*dble(iby))/dble(nby) 1249 2 continue 1250 do 3 ibz=1,nbz 1251 ibownr(ibz,3)=((ibz-1)*npz)/nbz 1252 boxsiz(ibz,3)=(box(3)*dble(ibz))/dble(nbz) 1253 3 continue 1254 if(iand(idebug,1).eq.1) then 1255 write(lfndbg,8003) nbx,nby,nbz,maxbox,npx,npy,npz 1256 8003 format('ibownr in sp_decomp',7i5) 1257 write(lfndbg,8004) (ibownr(ibx,1),ibx=1,nbx) 1258 8004 format('ibownr x',/,(20i5)) 1259 write(lfndbg,8005) (ibownr(iby,2),iby=1,nby) 1260 8005 format('ibownr y',/,(20i5)) 1261 write(lfndbg,8006) (ibownr(ibz,3),ibz=1,nbz) 1262 8006 format('ibownr z',/,(20i5)) 1263 write(lfndbg,8002) 1264 8002 format('boxlist') 1265 endif 1266c 1267c determine neighboring nodes and store in neighb(27,2) 1268c such that: 1269c 1270c neighb(n,1) is the n-th neighbor of node me 1271c neighb(n,2) is the node of which node me is the n-th neighbor 1272c 1273c a value of -1 indicates that such node does not exist 1274c 1275 do 4 i=1,27 1276 neighb(i,1)=-1 1277 neighb(i,2)=-1 1278 4 continue 1279c 1280 do 5 ix=1,3 1281 jx=mex+ix-2 1282 if(npbtyp.gt.0) then 1283 if(npx.gt.2.and.jx.lt.0) jx=jx+npx 1284 if(npx.gt.2.and.jx.ge.npx) jx=jx-npx 1285 endif 1286 if(jx.ge.0.and.jx.lt.npx) then 1287 do 6 iy=1,3 1288 jy=mey+iy-2 1289 if(npbtyp.gt.0) then 1290 if(npy.gt.2.and.jy.lt.0) jy=jy+npy 1291 if(npy.gt.2.and.jy.ge.npy) jy=jy-npy 1292 endif 1293 if(jy.ge.0.and.jy.lt.npy) then 1294 do 7 iz=1,3 1295 jz=mez+iz-2 1296 if(npbtyp.eq.1) then 1297 if(npz.gt.2.and.jz.lt.0) jz=jz+npz 1298 if(npz.gt.2.and.jz.ge.npz) jz=jz-npz 1299 endif 1300 if(jz.ge.0.and.jz.lt.npz) then 1301 jnode=((jz*npy)+jy)*npx+jx 1302 neighb(3*(3*(ix-1)+(iy-1))+iz,1)=jnode 1303 neighb(3*(3*(3-ix)+(3-iy))+4-iz,2)=jnode 1304 endif 1305 7 continue 1306 endif 1307 6 continue 1308 endif 1309 5 continue 1310c 1311 nbiown=30 1312 do 8 ibx=1,nbx 1313 do 9 iby=1,nby 1314 do 10 ibz=1,nbz 1315 nrnod=(ibownr(ibz,3)*npy+ibownr(iby,2))*npx+ibownr(ibx,1) 1316 if(me.eq.nrnod) nbiown=nbiown+1 1317 10 continue 1318 9 continue 1319 8 continue 1320c 1321 do 11 j=1,27 1322 do 12 i=1,np 1323 iburen(i,j,1)=0 1324 iburen(i,j,2)=0 1325 12 continue 1326 11 continue 1327 do 13 i=1,np 1328 ibindx(i)=0 1329 13 continue 1330 do 14 j=1,27 1331 iburen(me+1,j,1)=neighb(j,1) 1332 iburen(me+1,j,2)=neighb(j,2) 1333 if(neighb(j,1).ge.0) ibindx(neighb(j,1)+1)=j 1334 14 continue 1335c 1336 call ga_igop(msp_27,iburen,np*54,'+') 1337c 1338c broadcast the maximum number of sub-boxes per node to all nodes 1339c 1340 if(nbiown.gt.mbox) call md_abort('Error in mbox',0) 1341c 1342 return 1343 end 1344 subroutine sp_nrnode 1345c 1346 implicit none 1347c 1348#include "sp_common.fh" 1349c 1350 integer ix,iy,iz,npt,i,j,k,l 1351c 1352c this routine distributes the available processes in the 1353c Cartesian directions 1354c 1355c npx : number of processes in x direction 1356c npy : number of processes in y direction 1357c npz : number of processes in z direction 1358c 1359c determine node dimensions 1360c 1361 if(npx*npy*npz.ne.np) then 1362 if(npx+npy+npz.gt.0) 1363 + call md_abort('Specified npx*npy*npz ne np',0) 1364c 1365 npt=0 1366 do 1 i=1,np 1367 do 2 j=i,np 1368 do 3 k=j,np 1369 l=i*j*k 1370 if(l.eq.np) then 1371 if(l.eq.npt) then 1372 if(k.gt.npz) goto 3 1373 if(i+j+k.lt.npx+npy+npz) then 1374 npx=i 1375 npy=j 1376 npz=k 1377 endif 1378 goto 3 1379 else 1380 npt=np 1381 npx=i 1382 npy=j 1383 npz=k 1384 endif 1385 endif 1386 3 continue 1387 2 continue 1388 1 continue 1389c 1390 if(npx*npy*npz.ne.np) call md_abort('nrnode: code error',0) 1391 endif 1392c 1393c determine processor location of me 1394c 1395 do 4 ix=1,npx 1396 do 5 iy=1,npy 1397 do 6 iz=1,npz 1398 if(me.eq.((iz-1)*npy+(iy-1))*npx+(ix-1)) then 1399 mex=ix-1 1400 mey=iy-1 1401 mez=iz-1 1402 endif 1403 6 continue 1404 5 continue 1405 4 continue 1406c 1407 return 1408 end 1409 subroutine sp_initip(ibownr,ipl) 1410c 1411 implicit none 1412c 1413#include "sp_common.fh" 1414#include "global.fh" 1415c 1416 integer ibownr(maxbox,3),ipl(mbox,mip2) 1417 integer i,j,ibx,iby,ibz,nrbox,nrnod,nbiown 1418 integer il,ih,jl,jh 1419c 1420 do 1 j=1,mip2 1421 do 2 i=1,30 1422 ipl(i,j)=0 1423 2 continue 1424 1 continue 1425c 1426 if(iand(idebug,1).eq.1) then 1427 write(lfndbg,8003) nbx,nby,nbz,maxbox 1428 8003 format('ibownr in sp_initip',4i5) 1429 write(lfndbg,8004) (ibownr(ibx,1),ibx=1,nbx) 1430 8004 format('ibownr x',/,(20i5)) 1431 write(lfndbg,8005) (ibownr(iby,2),iby=1,nby) 1432 8005 format('ibownr y',/,(20i5)) 1433 write(lfndbg,8006) (ibownr(ibz,3),ibz=1,nbz) 1434 8006 format('ibownr z',/,(20i5)) 1435 call util_flush(lfndbg) 1436 if(iand(idebug,2).eq.2) write(lfndbg,8002) 1437 8002 format('boxlist') 1438 endif 1439 nbiown=0 1440 do 3 ibx=1,nbx 1441 do 4 iby=1,nby 1442 do 5 ibz=1,nbz 1443 nrbox=((ibz-1)*nby+iby-1)*nbx+ibx-1 1444 nrnod=(ibownr(ibz,3)*npy+ibownr(iby,2))*npx+ibownr(ibx,1) 1445 if(iand(idebug,2).eq.2) then 1446 write(lfndbg,8001) ibx,iby,ibz,nrbox,nrnod 1447 8001 format(5i5) 1448 call util_flush(lfndbg) 1449 endif 1450 if(me.eq.nrnod) then 1451 nbiown=nbiown+1 1452 ipl(30+nbiown,1)=nrbox 1453 ipl(30+nbiown,2)=0 1454 ipl(30+nbiown,3)=0 1455 ipl(30+nbiown,4)=0 1456 ipl(30+nbiown,5)=0 1457 endif 1458 5 continue 1459 4 continue 1460 3 continue 1461 ipl(1,1)=nbiown 1462c 1463 call ga_distribution(ga_ip,me,il,ih,jl,jh) 1464 call ga_put(ga_ip,il,ih,jl,jh,ipl,mbox) 1465c 1466 if(iand(idebug,2).eq.2) then 1467 write(lfndbg,8000) (i,ipl(30+i,1),i=1,ipl(1,1)) 1468 8000 format('ipl',/,(2i5)) 1469 call util_flush(lfndbg) 1470 endif 1471c 1472 return 1473 end 1474 subroutine sp_update_i(numsa,isl,numwm,iwl) 1475c 1476 implicit none 1477c 1478#include "sp_common.fh" 1479#include "mafdecls.fh" 1480#include "global.fh" 1481#include "msgids.fh" 1482c 1483 integer numsa,isl(msa,mis2) 1484 integer numwm,iwl(mwm,miw2) 1485c 1486 call sp_upd_i(numsa,isl,int_mb(i_pack),numwm,iwl,int_mb(i_packw)) 1487c 1488 return 1489 end 1490 subroutine sp_upd_i(numsa,isl,islp,numwm,iwl,iwlp) 1491c 1492 implicit none 1493c 1494#include "sp_common.fh" 1495#include "global.fh" 1496#include "msgids.fh" 1497c 1498 integer numsa,isl(msa,mis2) 1499 integer numwm,iwl(mwm,miw2) 1500 integer iwlp(mwm,npackw),islp(msa,npack) 1501c 1502 integer il,ih,jl,jh 1503c 1504 if(numsa.gt.0) then 1505 call ga_distribution(ga_is,me,il,ih,jl,jh) 1506 if(npack.eq.0) then 1507 call ga_put(ga_is,il,il+numsa-1,jl,jh,isl,msa) 1508 else 1509 call sp_pack(numsa,isl,islp) 1510 call ga_put(ga_is,il,il+numsa-1,jl,jl+npack-1,islp,msa) 1511 endif 1512 endif 1513c 1514 if(numwm.gt.0) then 1515 call ga_distribution(ga_iw,me,il,ih,jl,jh) 1516 if(npackw.eq.0) then 1517 call ga_put(ga_iw,il,il+numwm-1,jl,jh,iwl,mwm) 1518 else 1519 call sp_packw(numwm,iwl,iwlp) 1520 call ga_put(ga_iw,il,il+numwm-1,jl,jl+npackw-1,iwlp,mwm) 1521 endif 1522 endif 1523c 1524 call ga_sync() 1525c 1526 return 1527 end 1528 subroutine sp_numbb(ibownr,boxsiz) 1529c 1530 implicit none 1531c 1532#include "sp_common.fh" 1533#include "mafdecls.fh" 1534#include "msgids.fh" 1535#include "global.fh" 1536#include "util.fh" 1537c 1538 integer ibownr(maxbox,3) 1539 real*8 boxsiz(maxbox,3) 1540 logical lside,leven 1541c 1542 integer ibx,iby,ibz,ipx,ipy,ipz,ibox,inode 1543 integer jbx,jby,jbz,kbx,kby,kbz,jpx,jpy,jpz 1544 integer ilx,ihx,ily,ihy,ilz,ihz 1545 integer jbox,jnode,mbblb,i 1546 real*8 dx,dxtmp,dy,dytmp,dz,dztmp,dist2 1547 character*255 string 1548c 1549 if(iand(idebug,2).eq.2) then 1550 write(lfndbg,'(a,i6)') 'boxsiz in sp_numbb' 1551 write(lfndbg,'(6f12.6)') (boxsiz(i,1),i=1,nbx) 1552 write(lfndbg,'(6f12.6)') (boxsiz(i,2),i=1,nby) 1553 write(lfndbg,'(6f12.6)') (boxsiz(i,3),i=1,nbz) 1554 call util_flush(lfndbg) 1555 endif 1556c 1557c determine size of box-box pairlist 1558c ---------------------------------- 1559c 1560 mbbl=0 1561 mbblb=0 1562 do 1 ibx=0,nbx-1 1563 ipx=ibownr(ibx+1,1) 1564 do 2 iby=0,nby-1 1565 ipy=ibownr(iby+1,2) 1566 do 3 ibz=0,nbz-1 1567 ipz=ibownr(ibz+1,3) 1568 ibox=(ibz*nby+iby)*nbx+ibx 1569 inode=(ipz*npy+ipy)*npx+ipx 1570 if(inode.eq.me) then 1571 do 4 jbx=0,nbx-1 1572 kbx=jbx-ibx 1573 jpx=ibownr(jbx+1,1) 1574c 1575 if(ibx.le.jbx) then 1576 ilx=ibx 1577 ihx=jbx 1578 else 1579 ilx=jbx 1580 ihx=ibx 1581 endif 1582c 1583 dx=zero 1584 if(ibx.ne.jbx) then 1585 dx=boxsiz(ihx,1)-boxsiz(ilx+1,1) 1586 if(npbtyp.gt.0) then 1587 dxtmp=zero 1588 if(ilx.gt.0) dxtmp=boxsiz(ilx,1) 1589 if(ihx.lt.nbx-1) dxtmp=dxtmp-boxsiz(ihx+1,1)+box(1) 1590 if(dxtmp.lt.dx) dx=dxtmp 1591 if(kbx.gt.0.and.kbx.gt.iabs(kbx-nbx)) kbx=kbx-nbx 1592 if(kbx.lt.0.and.-kbx.gt.iabs(kbx+nbx)) kbx=kbx+nbx 1593 endif 1594 endif 1595c 1596 do 5 jby=0,nby-1 1597 kby=jby-iby 1598 jpy=ibownr(jby+1,2) 1599c 1600 if(iby.le.jby) then 1601 ily=iby 1602 ihy=jby 1603 else 1604 ily=jby 1605 ihy=iby 1606 endif 1607c 1608 dy=zero 1609 if(iby.ne.jby) then 1610 dy=boxsiz(ihy,2)-boxsiz(ily+1,2) 1611 if(npbtyp.gt.0) then 1612 dytmp=zero 1613 if(ily.gt.0) dytmp=boxsiz(ily,2) 1614 if(ihy.lt.nby-1) dytmp=dytmp-boxsiz(ihy+1,2)+box(2) 1615 if(dytmp.lt.dy) dy=dytmp 1616 if(kby.gt.0.and.kby.gt.iabs(kby-nby)) kby=kby-nby 1617 if(kby.lt.0.and.-kby.gt.iabs(kby+nby)) kby=kby+nby 1618 endif 1619 endif 1620c 1621 do 6 jbz=0,nbz-1 1622 kbz=jbz-ibz 1623 jpz=ibownr(jbz+1,3) 1624c 1625 if(ibz.le.jbz) then 1626 ilz=ibz 1627 ihz=jbz 1628 else 1629 ilz=jbz 1630 ihz=ibz 1631 endif 1632c 1633 dz=zero 1634 if(ibz.ne.jbz) then 1635 dz=boxsiz(ihz,3)-boxsiz(ilz+1,3) 1636 if(npbtyp.eq.1) then 1637 dztmp=zero 1638 if(ilz.gt.0) dztmp=boxsiz(ilz,3) 1639 if(ihz.lt.nbz-1) dztmp=dztmp-boxsiz(ihz+1,3)+box(3) 1640 if(dztmp.lt.dz) dz=dztmp 1641 if(kbz.gt.0.and.kbz.gt.iabs(kbz-nbz)) kbz=kbz-nbz 1642 if(kbz.lt.0.and.-kbz.gt.iabs(kbz+nbz)) kbz=kbz+nbz 1643 endif 1644 endif 1645c 1646 jbox=(jbz*nby+jby)*nbx+jbx 1647 jnode=(jpz*npy+jpy)*npx+jpx 1648c 1649c determine orientation jbox in relation to ibox 1650c 1651c lside is true if 1652c 1653c i: 0 j: 0 k: + 1654c i: 0 j: + k:-0+ 1655c i: + j:-0+ k:-0+ 1656c 1657 lside=(kbx.eq.0.and.kby.eq.0.and.kbz.ge.0) 1658 + .or.(kbx.eq.0.and.kby.gt.0) .or. kbx.gt.0 1659c 1660c determine if ibox is identical to jbox 1661c 1662c lsame=kbx.eq.0.and.kby.eq.0.and.kbz.eq.0 1663c 1664c determine if difference in box numbers is even or odd 1665c 1666 leven=2*(iabs(ibox-jbox)/2).eq.iabs(ibox-jbox) 1667c 1668c calculate the distance between the two boxes 1669c 1670 if(nbxtyp.eq.1) then 1671 dist2= 1672 + (vlat(1,1)*dx/box(1)+vlat(1,2)*dy/box(2)+vlat(1,3)*dz/box(3))**2+ 1673 + (vlat(2,1)*dx/box(1)+vlat(2,2)*dy/box(2)+vlat(2,3)*dz/box(3))**2+ 1674 + (vlat(3,1)*dx/box(1)+vlat(3,2)*dy/box(2)+vlat(3,3)*dz/box(3))**2 1675 else 1676 dist2=dx*dx+dy*dy+dz*dz 1677 endif 1678c 1679c keep half of the box pairs 1680c 1681c this test also appears in sp_numbb 1682c any changes need to be made in both routines 1683c 1684 if((inode.eq.jnode.and.ibox.ge.jbox).or.(inode.ne.jnode.and. 1685 + ((lside.and.leven).or.(.not.lside.and..not.leven)))) then 1686c 1687c keep only those within maximum cutoff distance 1688c 1689 if(rbbl*rbbl.gt.dist2) mbbl=mbbl+1 1690 endif 1691 if(rbbl*rbbl.gt.dist2) mbblb=mbblb+1 1692 6 continue 1693 5 continue 1694 4 continue 1695 endif 1696 3 continue 1697 2 continue 1698 1 continue 1699c 1700 mbbl=mbbl+1 1701 if(mbblb+1.gt.mbbl) mbbl=mbblb+1 1702c 1703 if(np.gt.1) call ga_igop(msp_20,mbbl,1,'max') 1704 if(me.eq.0) then 1705 if(util_print('distribution',print_debug)) then 1706 write(lfnout,2001) mbbl 1707 2001 format(' Dimension of the cell-cell list is ',i7) 1708 endif 1709 endif 1710c 1711c determine size of box-box pairlist assuming minimum cell sizes 1712c -------------------------------------------------------------- 1713c 1714 mbbl=0 1715 mbblb=0 1716 do 11 ibx=0,nbx-1 1717 ipx=ibownr(ibx+1,1) 1718 do 12 iby=0,nby-1 1719 ipy=ibownr(iby+1,2) 1720 do 13 ibz=0,nbz-1 1721 ipz=ibownr(ibz+1,3) 1722 ibox=(ibz*nby+iby)*nbx+ibx 1723 inode=(ipz*npy+ipy)*npx+ipx 1724 if(inode.eq.me) then 1725 do 14 jbx=0,nbx-1 1726 kbx=jbx-ibx 1727 jpx=ibownr(jbx+1,1) 1728c 1729 if(ibx.le.jbx) then 1730 ilx=ibx 1731 ihx=jbx 1732 else 1733 ilx=jbx 1734 ihx=ibx 1735 endif 1736c 1737 dx=zero 1738 if(ibx.ne.jbx) then 1739 dx=dble(ihx-ilx-1)*bxmin 1740 if(npbtyp.gt.0) then 1741 dxtmp=zero 1742 if(ilx.gt.0) dxtmp=dble(ilx)*bxmin 1743 if(ihx.lt.nbx-1) dxtmp=dxtmp+dble(nbx-ihx-1)*bxmin 1744 if(dxtmp.lt.dx) dx=dxtmp 1745 if(kbx.gt.0.and.kbx.gt.iabs(kbx-nbx)) kbx=kbx-nbx 1746 if(kbx.lt.0.and.-kbx.gt.iabs(kbx+nbx)) kbx=kbx+nbx 1747 endif 1748 endif 1749c 1750 do 15 jby=0,nby-1 1751 kby=jby-iby 1752 jpy=ibownr(jby+1,2) 1753c 1754 if(iby.le.jby) then 1755 ily=iby 1756 ihy=jby 1757 else 1758 ily=jby 1759 ihy=iby 1760 endif 1761c 1762 dy=zero 1763 if(iby.ne.jby) then 1764 dy=dble(ihy-ily-1)*bymin 1765 if(npbtyp.gt.0) then 1766 dytmp=zero 1767 if(ily.gt.0) dytmp=dble(ily)*bymin 1768 if(ihy.lt.nby-1) dytmp=dytmp+dble(nby-ihy-1)*bymin 1769 if(dytmp.lt.dy) dy=dytmp 1770 if(kby.gt.0.and.kby.gt.iabs(kby-nby)) kby=kby-nby 1771 if(kby.lt.0.and.-kby.gt.iabs(kby+nby)) kby=kby+nby 1772 endif 1773 endif 1774c 1775 do 16 jbz=0,nbz-1 1776 kbz=jbz-ibz 1777 jpz=ibownr(jbz+1,3) 1778c 1779 if(ibz.le.jbz) then 1780 ilz=ibz 1781 ihz=jbz 1782 else 1783 ilz=jbz 1784 ihz=ibz 1785 endif 1786c 1787 dz=zero 1788 if(ibz.ne.jbz) then 1789 dz=dble(ihz-ilz-1)*bzmin 1790 if(npbtyp.eq.1) then 1791 dztmp=zero 1792 if(ilz.gt.0) dztmp=dble(ilz)*bzmin 1793 if(ihz.lt.nbz-1) dztmp=dztmp-dble(nbz-ihz-1)*bzmin 1794 if(dztmp.lt.dz) dz=dztmp 1795 if(kbz.gt.0.and.kbz.gt.iabs(kbz-nbz)) kbz=kbz-nbz 1796 if(kbz.lt.0.and.-kbz.gt.iabs(kbz+nbz)) kbz=kbz+nbz 1797 endif 1798 endif 1799c 1800 jbox=(jbz*nby+jby)*nbx+jbx 1801 jnode=(jpz*npy+jpy)*npx+jpx 1802c 1803c determine orientation jbox in relation to ibox 1804c 1805c lside is true if 1806c 1807c i: 0 j: 0 k: + 1808c i: 0 j: + k:-0+ 1809c i: + j:-0+ k:-0+ 1810c 1811 lside=(kbx.eq.0.and.kby.eq.0.and.kbz.ge.0) 1812 + .or.(kbx.eq.0.and.kby.gt.0) .or. kbx.gt.0 1813c 1814c determine if ibox is identical to jbox 1815c 1816c lsame=kbx.eq.0.and.kby.eq.0.and.kbz.eq.0 1817c 1818c determine if difference in box numbers is even or odd 1819c 1820 leven=2*(iabs(ibox-jbox)/2).eq.iabs(ibox-jbox) 1821c 1822c calculate the distance between the two boxes 1823c 1824 if(nbxtyp.eq.1) then 1825 dist2= 1826 + (vlat(1,1)*dx/box(1)+vlat(1,2)*dy/box(2)+vlat(1,3)*dz/box(3))**2+ 1827 + (vlat(2,1)*dx/box(1)+vlat(2,2)*dy/box(2)+vlat(2,3)*dz/box(3))**2+ 1828 + (vlat(3,1)*dx/box(1)+vlat(3,2)*dy/box(2)+vlat(3,3)*dz/box(3))**2 1829 else 1830 dist2=dx*dx+dy*dy+dz*dz 1831 endif 1832c 1833c keep half of the box pairs 1834c 1835c this test also appears in sp_numbb 1836c any changes need to be made in both routines 1837c 1838 if((inode.eq.jnode.and.ibox.ge.jbox).or.(inode.ne.jnode.and. 1839 + ((lside.and.leven).or.(.not.lside.and..not.leven)))) then 1840c 1841c keep only those within maximum cutoff distance 1842c 1843 if(rbbl*rbbl.gt.dist2) mbbl=mbbl+1 1844 endif 1845 if(rbbl*rbbl.gt.dist2) mbblb=mbblb+1 1846 16 continue 1847 15 continue 1848 14 continue 1849 endif 1850 13 continue 1851 12 continue 1852 11 continue 1853c 1854 mbbl=mbbl+1 1855 if(mbblb+1.gt.mbbl) mbbl=mbblb+1 1856c 1857c 1858c allocate memory for the box-box list 1859c 1860 if(np.gt.1) call ga_igop(msp_20,mbbl,1,'max') 1861c 1862 if(me.eq.0) then 1863 if(util_print('distribution',print_debug)) then 1864 write(lfnout,2002) mbbl 1865 2002 format(' Dimension of the cell-cell list is ',i7) 1866 endif 1867 endif 1868c 1869 if(mbblp.eq.0) then 1870 mbbl=max(mbbl,mbbreq) 1871 if(.not.ma_push_get(mt_int,mbb2*mbbl,'bb',l_bb,i_bb)) 1872 + call md_abort('Failed to allocate memory for bb',0) 1873 mbblp=mbbl 1874 else 1875 if(mbbl.gt.mbblp) then 1876 write(string,1111) mbblp,mbbl,mbbreq 1877 1111 format('error: lbbl increase from ',i6,' to ',i6,'(',i6,')') 1878 call md_abort(string,me) 1879c call md_abort('lbbl increased beyond allocated memory',me) 1880 endif 1881 endif 1882c 1883 return 1884 end 1885 subroutine sp_listbb(ibownr,boxsiz,lbbl) 1886c 1887 implicit none 1888c 1889#include "sp_common.fh" 1890#include "global.fh" 1891#include "msgids.fh" 1892c 1893 integer ibownr(maxbox,3),lbbl(mbbl,mbb2) 1894 real*8 boxsiz(maxbox,3) 1895 logical lside,leven 1896c 1897 integer ibx,iby,ibz,ipx,ipy,ipz,ibox,inode 1898 integer jbx,jby,jbz,kbx,kby,kbz,jpx,jpy,jpz 1899 integer ilx,ihx,ily,ihy,ilz,ihz 1900 integer i,j,jbox,jnode,ltemp 1901 real*8 dx,dxtmp,dy,dytmp,dz,dztmp,dist2 1902c 1903c Construction of the box-box pairlist 1904c 1905 nbbl=0 1906 do 1 ibx=0,nbx-1 1907 ipx=ibownr(ibx+1,1) 1908 do 2 iby=0,nby-1 1909 ipy=ibownr(iby+1,2) 1910 do 3 ibz=0,nbz-1 1911 ipz=ibownr(ibz+1,3) 1912 ibox=(ibz*nby+iby)*nbx+ibx 1913 inode=(ipz*npy+ipy)*npx+ipx 1914 if(inode.eq.me) then 1915 do 4 jbx=0,nbx-1 1916 kbx=jbx-ibx 1917 jpx=ibownr(jbx+1,1) 1918c 1919 if(ibx.le.jbx) then 1920 ilx=ibx 1921 ihx=jbx 1922 else 1923 ilx=jbx 1924 ihx=ibx 1925 endif 1926c 1927 dx=zero 1928 if(ibx.ne.jbx) then 1929 dx=boxsiz(ihx,1)-boxsiz(ilx+1,1) 1930 if(npbtyp.gt.0) then 1931 dxtmp=zero 1932 if(ilx.gt.0) dxtmp=boxsiz(ilx,1) 1933 if(ihx.lt.nbx-1) dxtmp=dxtmp-boxsiz(ihx+1,1)+box(1) 1934 if(dxtmp.lt.dx) dx=dxtmp 1935 if(kbx.gt.0.and.kbx.gt.iabs(kbx-nbx)) kbx=kbx-nbx 1936 if(kbx.lt.0.and.-kbx.gt.iabs(kbx+nbx)) kbx=kbx+nbx 1937 endif 1938 endif 1939c 1940 do 5 jby=0,nby-1 1941 kby=jby-iby 1942 jpy=ibownr(jby+1,2) 1943c 1944 if(iby.le.jby) then 1945 ily=iby 1946 ihy=jby 1947 else 1948 ily=jby 1949 ihy=iby 1950 endif 1951c 1952 dy=zero 1953 if(iby.ne.jby) then 1954 dy=boxsiz(ihy,2)-boxsiz(ily+1,2) 1955 if(npbtyp.gt.0) then 1956 dytmp=zero 1957 if(ily.gt.0) dytmp=boxsiz(ily,2) 1958 if(ihy.lt.nby-1) dytmp=dytmp-boxsiz(ihy+1,2)+box(2) 1959 if(dytmp.lt.dy) dy=dytmp 1960 if(kby.gt.0.and.kby.gt.iabs(kby-nby)) kby=kby-nby 1961 if(kby.lt.0.and.-kby.gt.iabs(kby+nby)) kby=kby+nby 1962 endif 1963 endif 1964c 1965 do 6 jbz=0,nbz-1 1966 kbz=jbz-ibz 1967 jpz=ibownr(jbz+1,3) 1968c 1969 if(ibz.le.jbz) then 1970 ilz=ibz 1971 ihz=jbz 1972 else 1973 ilz=jbz 1974 ihz=ibz 1975 endif 1976c 1977 dz=zero 1978 if(ibz.ne.jbz) then 1979 dz=boxsiz(ihz,3)-boxsiz(ilz+1,3) 1980 if(npbtyp.eq.1) then 1981 dztmp=zero 1982 if(ilz.gt.0) dztmp=boxsiz(ilz,3) 1983 if(ihz.lt.nbz-1) dztmp=dztmp-boxsiz(ihz+1,3)+box(3) 1984 if(dztmp.lt.dz) dz=dztmp 1985 if(kbz.gt.0.and.kbz.gt.iabs(kbz-nbz)) kbz=kbz-nbz 1986 if(kbz.lt.0.and.-kbz.gt.iabs(kbz+nbz)) kbz=kbz+nbz 1987 endif 1988 endif 1989c 1990 jbox=(jbz*nby+jby)*nbx+jbx 1991 jnode=(jpz*npy+jpy)*npx+jpx 1992c 1993c determine orientation jbox in relation to ibox 1994c 1995c lside is true if 1996c 1997c i: 0 j: 0 k: + 1998c i: 0 j: + k:-0+ 1999c i: + j:-0+ k:-0+ 2000c 2001 lside=(kbx.eq.0.and.kby.eq.0.and.kbz.ge.0) 2002 + .or.(kbx.eq.0.and.kby.gt.0) .or. kbx.gt.0 2003c 2004c determine if ibox is identical to jbox 2005c 2006c lsame=kbx.eq.0.and.kby.eq.0.and.kbz.eq.0 2007c 2008c determine if difference in box numbers is even or odd 2009c 2010 leven=2*(iabs(ibox-jbox)/2).eq.iabs(ibox-jbox) 2011c 2012c calculate the distance between the two boxes 2013c 2014 if(nbxtyp.eq.1) then 2015 dist2= 2016 + (vlat(1,1)*dx/box(1)+vlat(1,2)*dy/box(2)+vlat(1,3)*dz/box(3))**2+ 2017 + (vlat(2,1)*dx/box(1)+vlat(2,2)*dy/box(2)+vlat(2,3)*dz/box(3))**2+ 2018 + (vlat(3,1)*dx/box(1)+vlat(3,2)*dy/box(2)+vlat(3,3)*dz/box(3))**2 2019 else 2020 dist2=dx*dx+dy*dy+dz*dz 2021 endif 2022c 2023c keep half of the box pairs 2024c 2025c this test also appears in sp_numbb 2026c any changes need to be made in both routines 2027c 2028 if((inode.eq.jnode.and.ibox.ge.jbox).or. (inode.ne.jnode.and. 2029 + ((lside.and.leven).or.(.not.lside.and..not.leven)))) then 2030c 2031c keep only those within maximum cutoff distance 2032c 2033 if(rbbl*rbbl.gt.dist2) then 2034 nbbl=nbbl+1 2035 if(nbbl.gt.mbbl) call md_abort('Box-box list too small',mbbl) 2036 lbbl(nbbl,1)=jnode 2037 lbbl(nbbl,2)=jbox 2038 lbbl(nbbl,3)=ibox 2039 lbbl(nbbl,4)=0 2040 endif 2041 endif 2042 6 continue 2043 5 continue 2044 4 continue 2045 endif 2046 3 continue 2047 2 continue 2048 1 continue 2049 npprev=0 2050c 2051 nbbloc=0 2052 do 7 i=1,nbbl-1 2053 do 8 j=i+1,nbbl 2054 if((lbbl(i,1).ne.me.and.lbbl(j,1).eq.me).or. 2055 + (lbbl(i,1).gt.lbbl(j,1).and.lbbl(i,1).ne.me).or. 2056 + (lbbl(i,1).eq.lbbl(j,1).and.lbbl(i,2).gt.lbbl(j,2)).or. 2057 + (lbbl(i,1).eq.lbbl(j,1).and.lbbl(i,2).eq.lbbl(j,2).and. 2058 + lbbl(i,3).gt.lbbl(j,3))) then 2059 ltemp=lbbl(i,1) 2060 lbbl(i,1)=lbbl(j,1) 2061 lbbl(j,1)=ltemp 2062 ltemp=lbbl(i,2) 2063 lbbl(i,2)=lbbl(j,2) 2064 lbbl(j,2)=ltemp 2065 ltemp=lbbl(i,3) 2066 lbbl(i,3)=lbbl(j,3) 2067 lbbl(j,3)=ltemp 2068 endif 2069 8 continue 2070 if(lbbl(i,1).eq.me) nbbloc=i 2071 7 continue 2072 if(lbbl(nbbl,1).eq.me) nbbloc=nbbl 2073c 2074 nrempr=0 2075 if(nbget.ne.0.and.np.gt.1) then 2076 nrempr=1 2077 do 9 i=2,nbbl 2078 if(lbbl(i,2).ne.lbbl(i-1,2)) nrempr=nrempr+1 2079 9 continue 2080 call ga_igop(msp_29,nrempr,1,'max') 2081 nrempr=min(2*nrempr,nrempr+25) 2082 endif 2083c 2084 if(iand(idebug,2).eq.2) then 2085 write(lfndbg,8000) (i,(lbbl(i,j),j=1,3),i=1,nbbl) 2086 8000 format('lbbl',/,(4i5)) 2087 call util_flush(lfndbg) 2088 endif 2089c 2090 return 2091 end 2092 subroutine sp_finish() 2093c 2094 implicit none 2095c 2096#include "sp_common.fh" 2097#include "mafdecls.fh" 2098c 2099 if(.not.ma_pop_stack(l_xscr)) 2100 + call md_abort('Failed to deallocate xscr',0) 2101 if(.not.ma_pop_stack(l_bb)) 2102 + call md_abort('Failed to deallocate memory for bb',0) 2103c 2104 call sp_free() 2105c 2106 return 2107 end 2108 subroutine sp_initf(fw,fs,llng,iwz,isz,lpair) 2109c 2110 implicit none 2111c 2112#include "sp_common.fh" 2113#include "mafdecls.fh" 2114#include "global.fh" 2115c 2116 real*8 fw(mwm,3,mwa,2),fs(msa,3,2) 2117 integer iwz(mwm),isz(msa) 2118 logical llng,lpair 2119c 2120 integer i,j,k,l,m,il,ih,jl,jh 2121c 2122 llong=llng 2123c 2124 m=1 2125 if(llong) m=2 2126c 2127 do 1 l=1,m 2128 if(nwm.gt.0) then 2129 do 2 k=1,mwa 2130 do 3 j=1,3 2131 do 4 i=1,mwm 2132 fw(i,j,k,l)=zero 2133 4 continue 2134 3 continue 2135 2 continue 2136 endif 2137 if(nsa.gt.0) then 2138 do 5 j=1,3 2139 do 6 i=1,msa 2140 fs(i,j,l)=zero 2141 6 continue 2142 5 continue 2143 endif 2144 1 continue 2145c 2146 if(nwm.gt.0) then 2147 call ga_distribution(ga_w,me,il,ih,jl,jh) 2148 call ga_put(ga_w,il,ih,jl+6*mwa+3,jl+9*mwa+2,fw,mwm) 2149 if(llong) call ga_put(ga_w,il,ih,jl+9*mwa+3,jl+12*mwa+2, 2150 + fw(1,1,1,2),mwm) 2151 endif 2152 if(nsa.gt.0) then 2153 call ga_distribution(ga_s,me,il,ih,jl,jh) 2154 call ga_put(ga_s,il,ih,jl+6,jl+8,fs,msa) 2155 if(llong) call ga_put(ga_s,il,ih,jl+9,jl+11,fs(1,1,2),msa) 2156 endif 2157c 2158 if(lpair) then 2159 do 7 i=1,mwm 2160 iwz(i)=0 2161 7 continue 2162 do 8 i=1,msa 2163 isz(i)=0 2164 8 continue 2165 call ga_zero(ga_iwz) 2166 call ga_zero(ga_isz) 2167 endif 2168c 2169 return 2170 end 2171 subroutine sp_copyg(fw,fs) 2172c 2173 implicit none 2174c 2175#include "sp_common.fh" 2176#include "global.fh" 2177c 2178 real*8 fw(mwm,3,mwa),fs(msa,3) 2179c 2180 integer il,ih,jl,jh 2181c 2182 if(nwm.gt.0) then 2183 call ga_distribution(ga_w,me,il,ih,jl,jh) 2184 call ga_put(ga_w,il,ih,jl+6*mwa+3,jl+9*mwa+2,fw,mwm) 2185 endif 2186 if(nsa.gt.0) then 2187 call ga_distribution(ga_s,me,il,ih,jl,jh) 2188 call ga_put(ga_s,il,ih,jl+6,jl+8,fs,msa) 2189 endif 2190c 2191 return 2192 end 2193 subroutine sp_final(fw,fs,lpair,iwz,isz) 2194c 2195 implicit none 2196c 2197#include "sp_common.fh" 2198#include "mafdecls.fh" 2199#include "global.fh" 2200c 2201 real*8 fw(mwm,3,mwa,2),fs(msa,3,2) 2202 logical lpair 2203 integer iwz(mwm),isz(msa) 2204c 2205 integer i,il,ih,jl,jh 2206c 2207 if(np.gt.0) then 2208 if(nwm.gt.0) then 2209 call ga_distribution(ga_w,me,il,ih,jl,jh) 2210 call ga_acc(ga_w,il,ih,jl+6*mwa+3,jl+9*mwa+2,fw,mwm,one) 2211 if(llong) call ga_acc(ga_w,il,ih,jl+9*mwa+3,jl+12*mwa+2, 2212 + fw(1,1,1,2),mwm,one) 2213 call ga_get(ga_w,il,ih,jl+6*mwa+3,jl+9*mwa+2,fw,mwm) 2214 if(ltwin) call ga_get(ga_w,il,ih,jl+9*mwa+3,jl+12*mwa+2, 2215 + fw(1,1,1,2),mwm) 2216 endif 2217 if(nsa.gt.0) then 2218 call ga_distribution(ga_s,me,il,ih,jl,jh) 2219 call ga_acc(ga_s,il,ih,jl+6,jl+8,fs,msa,one) 2220 if(llong) call ga_acc(ga_s,il,ih,jl+9,jl+11,fs(1,1,2),msa,one) 2221 call ga_get(ga_s,il,ih,jl+6,jl+8,fs,msa) 2222 if(ltwin) call ga_get(ga_s,il,ih,jl+9,jl+11,fs(1,1,2),msa) 2223 endif 2224 endif 2225c 2226 if(lpair) then 2227 if(nwm.gt.0) then 2228 call ga_distribution(ga_iwz,me,il,ih,jl,jh) 2229 call ga_acc(ga_iwz,il,ih,1,1,iwz,mwm,1) 2230 call ga_get(ga_iwz,il,ih,1,1,iwz,mwm) 2231 do 1 i=1,nwmloc 2232 iwz(i)=min(1,iwz(i)) 2233 1 continue 2234 endif 2235 if(nsa.gt.0) then 2236 call ga_distribution(ga_isz,me,il,ih,jl,jh) 2237 call ga_acc(ga_isz,il,ih,1,1,isz,msa,1) 2238 call ga_get(ga_isz,il,ih,1,1,isz,msa) 2239 do 2 i=1,nsaloc 2240 isz(i)=min(1,isz(i)) 2241 2 continue 2242 endif 2243 endif 2244c 2245 return 2246 end 2247 subroutine sp_wrtrst(lfnrst,filrst,lveloc, 2248 + pres,temp,tempw,temps,iwl,xw,vw,fw,xwcr,isl,xs,vs,fs,xscr,prjct, 2249 + lseq) 2250c 2251 implicit none 2252c 2253#include "sp_common.fh" 2254#include "mafdecls.fh" 2255#include "global.fh" 2256c 2257 integer lfnrst 2258 character*255 filrst 2259 logical lveloc 2260 integer iwl(mwm,miw2),isl(msa,mis2),lseq(mseq) 2261 real*8 pres,temp,tempw,temps 2262 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3) 2263 real*8 xs(msa,3),vs(msa,3),fs(msa,3),xscr(msm,3) 2264 character*80 prjct 2265c 2266 integer lenscr 2267c 2268 project=prjct 2269c 2270 lenscr=ma_inquire_avail(mt_byte)/ 2271 + ((9*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+ 2272 + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1 2273 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx)) 2274 + call md_abort('Failed to allocate bx',0) 2275 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv)) 2276 + call md_abort('Failed to allocate bv',0) 2277 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bf',l_bf,i_bf)) 2278 + call md_abort('Failed to allocate bf',0) 2279 if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br)) 2280 + call md_abort('Failed to allocate br',0) 2281 if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi)) 2282 + call md_abort('Failed to allocate bi',0) 2283 if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n)) 2284 + call md_abort('Failed to allocate n',0) 2285c 2286 call sp_wtrst(lfnrst,filrst,lveloc,pres,temp,tempw,temps, 2287 + iwl,int_mb(i_packw),xw,vw,fw,xwcr,isl,int_mb(i_pack),xs,vs,fs, 2288 + xscr,int_mb(i_ipl),lenscr,int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv), 2289 + dbl_mb(i_bf),dbl_mb(i_br),int_mb(i_bi),dbl_mb(i_bx), 2290 + dbl_mb(i_bv),dbl_mb(i_bf),lseq) 2291c 2292 if(.not.ma_pop_stack(l_n)) 2293 + call md_abort('Failed to deallocate n',0) 2294 if(.not.ma_pop_stack(l_bi)) 2295 + call md_abort('Failed to deallocate bi',0) 2296 if(.not.ma_pop_stack(l_br)) 2297 + call md_abort('Failed to deallocate br',0) 2298 if(.not.ma_pop_stack(l_bf)) 2299 + call md_abort('Failed to deallocate bf',0) 2300 if(.not.ma_pop_stack(l_bv)) 2301 + call md_abort('Failed to deallocate bv',0) 2302 if(.not.ma_pop_stack(l_bx)) 2303 + call md_abort('Failed to deallocate bx',0) 2304c 2305 return 2306 end 2307 subroutine sp_wtrst(lfnrst,filrst,lveloc,pres,temp,tempw,temps, 2308 + iwl,iwlp,xw,vw,fw,xwcr,isl,islp,xs,vs,fs,xscr, 2309 + ipl,nb,ibw,bxw,bvw,bfw,brw,ibs,bxs,bvs,bfs,lseq) 2310c 2311 implicit none 2312c 2313#include "sp_common.fh" 2314#include "mafdecls.fh" 2315#include "global.fh" 2316c 2317 integer lfnrst,nb 2318 character*255 filrst 2319 logical lveloc 2320 real*8 pres,temp,tempw,temps 2321 integer iwl(mwm,miw2),isl(msa,mis2),lseq(mseq) 2322 integer iwlp(mwm,npackw),islp(msa,npack) 2323 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3) 2324 real*8 xs(msa,3),vs(msa,3),fs(msa,3),xscr(msm,3) 2325 integer ipl(mbox,mip2),ibw(nb),ibs(nb,2) 2326 real*8 bxw(nb,3,mwa),bvw(nb,3,mwa),bfw(nb,3,mwa),brw(nb,3) 2327 real*8 bxs(nb,3),bvs(nb,3),bfs(nb,3) 2328c 2329 integer i,j,k,node,ncyc,icyc,numw,nums,number,nwmn,nsan 2330 integer ilp,ihp,jlp,jhp,ili,ihi,jli,jhi,ilw,ihw,jlw,jhw 2331 integer ils,ihs,jls,jhs 2332 character*10 rdate,rtime 2333 character*18 user 2334#ifdef USE_POSIXF 2335 integer ilen,ierror 2336#endif 2337 integer idyn,idynp,ihop 2338 logical lforces 2339c 2340 lforces=iguide.ne.0 2341c 2342 if(ga_nodeid().eq.0) then 2343c 2344 call swatch(rdate,rtime) 2345#ifdef USE_POSIXF 2346 call pxfgetlogin(user, ilen, ierror) 2347#else 2348 call getlog(user) 2349#endif 2350 if(user(18:18).ne.' ') user=' ' 2351c 2352 rewind(lfnrst) 2353 write(lfnrst,1000) 2354 1000 format('Restart file',/,' ',/,' ') 2355 write(lfnrst,1001) 4.2,rdate,rtime,nhist,lforces 2356 1001 format(f12.6,2a10,i5,4x,l1) 2357 hist(nhist)(1:18)=user 2358 hist(nhist)(19:28)=rdate 2359 hist(nhist)(29:48)=rtime 2360 hist(nhist)(49:108)=project(1:60) 2361 do 10 i=1,nhist 2362 write(lfnrst,1009) hist(i) 2363 1009 format(a) 2364 10 continue 2365 write(lfnrst,1002) npbtyp,nbxtyp,rsgm,((vlat(i,j),j=1,3),i=1,3) 2366 1002 format(2i5,f12.6,/,(3f12.6)) 2367 write(lfnrst,1003) pres 2368 1003 format(1pe12.5) 2369 write(lfnrst,1004) temp,tempw,temps 2370 1004 format(3f12.6) 2371 write(lfnrst,1005) nwm,nwa,nsm,nsa,nwmc,nsf,nseq,0,0 2372 1005 format(7i10,2i5) 2373c 2374 if(nwm.gt.0) then 2375 number=0 2376 ncyc=nwm/nb+1 2377 numw=nb 2378 do 1 icyc=1,ncyc 2379 if(nwm-number.lt.numw) numw=nwm-number 2380c 2381c begin test code 10/31/2001 2382c initialize ibw to check that all atoms have been received 2383c 2384 do 1112 i=1,nb 2385 ibw(i)=-1 2386 1112 continue 2387c 2388c end test code 2389c 2390 do 2 node=np-1,0,-1 2391 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 2392 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 2393 nwmn=ipl(1,2) 2394 if(nwmn.gt.0) then 2395 call ga_distribution(ga_iw,node,ili,ihi,jli,jhi) 2396 if(npackw.eq.0) then 2397 call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+lwdyn-1,iwl,mwm) 2398 else 2399 call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm) 2400 call sp_unpackw(nwmn,iwl,iwlp) 2401 endif 2402 call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw) 2403 call ga_get(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm) 2404 if(lveloc) 2405 + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm) 2406 if(lforces) 2407 + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa+3,jlw+9*mwa+2,fw,mwm) 2408 call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa,jlw+6*mwa+2,xwcr,mwm) 2409 do 3 i=1,nwmn 2410 j=iwl(i,lwgmn)-number 2411 if(j.gt.0.and.j.le.numw) then 2412 do 4 k=1,nwa 2413 bxw(j,1,k)=xw(i,1,k) 2414 bxw(j,2,k)=xw(i,2,k) 2415 bxw(j,3,k)=xw(i,3,k) 2416 bvw(j,1,k)=vw(i,1,k) 2417 bvw(j,2,k)=vw(i,2,k) 2418 bvw(j,3,k)=vw(i,3,k) 2419 if(lforces) then 2420 bfw(j,1,k)=fw(i,1,k) 2421 bfw(j,2,k)=fw(i,2,k) 2422 bfw(j,3,k)=fw(i,3,k) 2423 endif 2424 4 continue 2425 brw(j,1)=xwcr(i,1) 2426 brw(j,2)=xwcr(i,2) 2427 brw(j,3)=xwcr(i,3) 2428 ibw(j)=iwl(i,lwdyn) 2429 endif 2430 3 continue 2431 endif 2432 2 continue 2433 do 5 i=1,numw 2434 if(lveloc) then 2435 write(lfnrst,1006) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa) 2436 else 2437 write(lfnrst,1006) ((bxw(i,j,k),j=1,3),(zero,j=1,3),k=1,nwa) 2438 endif 2439 1006 format(2x,6f13.8) 2440 if(lforces) write(lfnrst,1106) ((bfw(i,j,k),j=1,3),k=1,nwa) 2441 1106 format(2x,6e13.6) 2442 idyn=iand(ibw(i),12)/4 2443 idynp=iand(ibw(i),3) 2444 write(lfnrst,1007) idynp,idyn,(brw(i,k),k=1,3) 2445 1007 format(2i1,3f13.8) 2446c 2447c begin test code 10/31/2001 2448c check if al atoms have been received 2449c 2450 if(ibw(i).lt.0) 2451 + call md_abort('Missing solvent in wtrst',i) 2452c 2453c end test code 2454c 2455 5 continue 2456 number=number+numw 2457 1 continue 2458 endif 2459c 2460 if(nsa.gt.0) then 2461 number=0 2462 ncyc=nsa/nb+1 2463 nums=nb 2464 do 6 icyc=1,ncyc 2465 if(nsa-number.lt.nums) nums=nsa-number 2466c 2467c begin test code 10/31/2001 2468c initialize ibw to check that all atoms have been received 2469c 2470 do 1117 i=1,nb 2471 ibs(i,1)=-1 2472 ibs(i,2)=0 2473 1117 continue 2474c 2475c end test code 2476c 2477 do 7 node=np-1,0,-1 2478 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 2479 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 2480 nsan=ipl(2,2) 2481 if(nsan.gt.0) then 2482 call ga_distribution(ga_is,node,ili,ihi,jli,jhi) 2483 if(npack.eq.0) then 2484 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+lsdyn-1,isl,msa) 2485 else 2486 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa) 2487 call sp_unpack(nsan,isl,islp) 2488 endif 2489 call ga_distribution(ga_s,node,ils,ihs,jls,jhs) 2490 call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa) 2491 if(lveloc) call ga_get(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa) 2492 if(lforces) call ga_get(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa) 2493 do 8 i=1,nsan 2494 j=isl(i,lsgan)-number 2495 if(j.gt.0.and.j.le.nums) then 2496 bxs(j,1)=xs(i,1) 2497 bxs(j,2)=xs(i,2) 2498 bxs(j,3)=xs(i,3) 2499 bvs(j,1)=vs(i,1) 2500 bvs(j,2)=vs(i,2) 2501 bvs(j,3)=vs(i,3) 2502 if(lforces) then 2503 bfs(j,1)=fs(i,1) 2504 bfs(j,2)=fs(i,2) 2505 bfs(j,3)=fs(i,3) 2506 endif 2507 ibs(j,1)=isl(i,lsdyn) 2508 ibs(j,2)=isl(i,lshop) 2509 endif 2510 8 continue 2511 endif 2512 7 continue 2513 do 9 i=1,nums 2514 idyn=iand(ibs(i,1),12)/4 2515 idynp=iand(ibs(i,1),3) 2516 ihop=ibs(i,2) 2517 if(iand(ihop,1).eq.1) then 2518 ihop=-(ihop/2) 2519 else 2520 ihop=ihop/2 2521 endif 2522 if(lveloc) then 2523 write(lfnrst,1008) idynp,idyn,(bxs(i,j),j=1,3),(bvs(i,j),j=1,3), 2524 + ihop 2525 else 2526 write(lfnrst,1008) idynp,idyn,(bxs(i,j),j=1,3),(zero,j=1,3),ihop 2527 endif 2528 1008 format(2i1,6f13.8,i5) 2529 if(lforces) write(lfnrst,1108) (bfs(i,j),j=1,3) 2530 1108 format(2x,3e13.6) 2531c 2532c begin test code 10/31/2001 2533c check if al atoms have been received 2534c 2535 if(ibs(i,1).lt.0) 2536 + call md_abort('Missing solute atom in wtrst',i) 2537c 2538c end test code 2539c 2540 9 continue 2541 number=number+nums 2542 6 continue 2543 endif 2544c 2545 if(nsm.gt.0) then 2546 do 21 i=1,nsm 2547 write(lfnrst,1109) (xscr(i,j),j=1,3) 2548 1109 format(2x,3f13.8) 2549 21 continue 2550 endif 2551c 2552 if(nseq.gt.0) then 2553 write(lfnrst,1013) (lseq(i),i=1,nseq) 2554 1013 format(20i3) 2555 endif 2556c 2557 endif 2558c 2559 return 2560 9999 continue 2561 call md_abort('Failed to open restart file',me) 2562 return 2563 end 2564 subroutine sp_print() 2565c 2566 implicit none 2567c 2568#include "sp_common.fh" 2569#include "mafdecls.fh" 2570c 2571 integer i_lcnt,l_lcnt 2572c 2573 if(me.eq.0) then 2574 write(lfnout,1000) 2575 1000 format(/,' DOMAIN DECOMPOSITION',/) 2576c 2577 write(lfnout,1001) np,npx,npy,npz 2578 1001 format(' Processor count ',i5,' =',i5,' x',i5,' x',i5) 2579 write(lfnout,1002) nbx*nby*nbz,nbx,nby,nbz 2580 1002 format(' Cell count ',i5,' =',i5,' x',i5,' x',i5) 2581c 2582 if(mod(nbx,npx)+mod(nby,npy)+mod(nbz,npz).ne.0) then 2583 write(lfnout,1003) 2584 1003 format(/,' WARNING: Inefficient distribution of cells over ', 2585 + 'processors') 2586 endif 2587c 2588 write(lfnout,1004) bxmin,bymin,bzmin 2589 1004 format(/,' Minimum cell size ',f12.6,2(' x',f12.6)) 2590c 2591 if(nred(1)+nred(2)+nred(3).gt.0) then 2592 write(lfnout,1005) 2593 1005 format(/,' Warning: ',/) 2594 if(nred(1).gt.0) write(lfnout,1006) 'x',nred(1) 2595 if(nred(2).gt.0) write(lfnout,1006) 'y',nred(2) 2596 if(nred(3).gt.0) write(lfnout,1006) 'z',nred(3) 2597 1006 format(' Reduced number of cells in ',a,'-dimension: ',i5) 2598 endif 2599c 2600 if(nable.eq.1) write(lfnout,1007) 2601 1007 format(/,' Read previous box pair list',/) 2602c 2603 if(nable.eq.2) write(lfnout,1008) 2604 1008 format(/,' Unable to read previous box pair list',/) 2605c 2606 endif 2607c 2608 if(.not.ma_push_get(mt_int,3*np,'lcnt',l_lcnt,i_lcnt)) 2609 + call md_abort('Failed to allocate memory for lcnt',0) 2610c 2611 call sp_prtcnt(int_mb(i_lcnt)) 2612c 2613 if(.not.ma_pop_stack(l_lcnt)) 2614 + call md_abort('Failed to deallocate lcnt',0) 2615c 2616 return 2617 end 2618 subroutine sp_prtcnt(lcnt) 2619c 2620 implicit none 2621c 2622#include "sp_common.fh" 2623#include "msgids.fh" 2624#include "global.fh" 2625c 2626 integer lcnt(np,3) 2627c 2628 integer i,j 2629c 2630 do 1 i=1,np 2631 lcnt(i,1)=0 2632 lcnt(i,2)=0 2633 lcnt(i,3)=0 2634 1 continue 2635c 2636 lcnt(me+1,1)=mbxloc 2637 lcnt(me+1,2)=nwmloc*nwa 2638 lcnt(me+1,3)=nsaloc 2639c 2640 if(np.gt.1) call ga_igop(msp_08,lcnt,3*np,'+') 2641c 2642 if(me.eq.0) then 2643 write(lfnout,1000) 2644 1000 format(/,' Initial distribution p:b(w+s)',/) 2645 write(lfnout,1001) (i-1,(lcnt(i,j),j=1,3),i=1,np) 2646 1001 format(4(3x,i4,':',i5,'(',i7,'+',i7,')')) 2647 write(lfnout,1002) mwm,msa 2648 1002 format(/,' Dimension workarrays solvent ',i6,/, 2649 + 22x,'solute ',i6) 2650 endif 2651c 2652 return 2653 end 2654 subroutine sp_printf(filtop,lfntop,isl,xs,fs, 2655 + npener,esa) 2656c 2657 implicit none 2658c 2659#include "sp_common.fh" 2660#include "mafdecls.fh" 2661#include "msgids.fh" 2662#include "global.fh" 2663c 2664 integer lfntop 2665 character*70 filtop 2666 integer isl(msa,mis2),npener 2667 real*8 xs(msa,3),fs(msa,3,2),esa(nsa,2) 2668c 2669 integer lenscr 2670c 2671 lenscr=ma_inquire_avail(mt_byte)/ 2672 + ((6*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+ 2673 + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1 2674 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx)) 2675 + call md_abort('Failed to allocate bx',0) 2676 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv)) 2677 + call md_abort('Failed to allocate bv',0) 2678 if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br)) 2679 + call md_abort('Failed to allocate br',0) 2680 if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi)) 2681 + call md_abort('Failed to allocate bi',0) 2682 if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n)) 2683 + call md_abort('Failed to allocate n',0) 2684c 2685 call sp_prt_s(filtop,lfntop, 2686 + int_mb(i_ipl),isl,int_mb(i_pack),xs,fs, 2687 + lenscr,int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv), 2688 + npener,esa) 2689c 2690 if(.not.ma_pop_stack(l_n)) 2691 + call md_abort('Failed to deallocate n',0) 2692 if(.not.ma_pop_stack(l_bi)) 2693 + call md_abort('Failed to deallocate bi',0) 2694 if(.not.ma_pop_stack(l_br)) 2695 + call md_abort('Failed to deallocate br',0) 2696 if(.not.ma_pop_stack(l_bv)) 2697 + call md_abort('Failed to deallocate bv',0) 2698 if(.not.ma_pop_stack(l_bx)) 2699 + call md_abort('Failed to deallocate bx',0) 2700c 2701 return 2702 end 2703 subroutine sp_prt_s(filtop,lfntop, 2704 + ipl,isl,islp,xs,fs,nb,ibs,bxs,bfs, 2705 + npener,esa) 2706c 2707 implicit none 2708c 2709#include "sp_common.fh" 2710#include "msgids.fh" 2711#include "global.fh" 2712c 2713 integer lfntop,nb,npener 2714 character*70 filtop 2715 integer ipl(mbox,mip2),isl(msa,mis2),islp(msa,npack),ibs(nb) 2716 real*8 xs(msa,3),fs(msa,3) 2717 real*8 bxs(nb,3),bfs(nb,3) 2718 real*8 esa(nsa,2) 2719c 2720 integer i,j,k,number,icyc,ncyc,node,nums,nsan,idyn,ism,iss 2721 integer ilp,ihp,jlp,jhp,ili,ihi,jli,jhi,ils,ihs,jls,jhs 2722 integer naw,nbw,nhw,ndw,now,ntw,nnw,nat,nqt,idum,npars 2723 character*1 cdum 2724 character*16 cat 2725c 2726 call ga_distribution(ga_ip,me,ilp,ihp,jlp,jhp) 2727 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 2728 nsan=ipl(2,2) 2729c 2730 if(nsan.gt.0) then 2731 call ga_distribution(ga_s,me,ils,ihs,jls,jhs) 2732 call ga_put(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa) 2733 endif 2734c 2735 call ga_sync() 2736c 2737 if(me.ne.0) return 2738c 2739 if(npener.le.0) then 2740 write(lfnout,1007) 2741 1007 format(//,' Solute coordinates and forces',//, 2742 + ' mol segment atom dt x y z ', 2743 + ' fx fy fz',/) 2744 else 2745 write(lfnout,1008) 2746 1008 format(//,' Solute coordinates, forces and energies',//, 2747 + ' mol segment atom dt x y z ', 2748 + ' fx fy fz enb', 2749 + ' eb ep',/) 2750 endif 2751c 2752 open(unit=lfntop,file=filtop(1:index(filtop,' ')-1), 2753 + status='old',form='formatted') 2754 rewind(lfntop) 2755 do 101 i=1,4 2756 read(lfntop,2001) cdum 2757 2001 format(a1) 2758 101 continue 2759 read(lfntop,2002) npars 2760 read(lfntop,2002) nat 2761 read(lfntop,2002) nqt 2762 read(lfntop,2002) nseq 2763 2002 format(i5) 2764 read(lfntop,2001) cdum 2765 do 1102 k=1,npars 2766 do 102 i=1,nat 2767 read(lfntop,2001) cdum 2768 do 103 j=i,nat 2769 read(lfntop,2001) cdum 2770 103 continue 2771 102 continue 2772 1102 continue 2773 do 104 i=1,nqt*npars 2774 read(lfntop,2001) cdum 2775 104 continue 2776 do 4104 i=1,nseq 2777 read(lfntop,2001) cdum 2778 4104 continue 2779 read(lfntop,2003) naw,nbw,nhw,ndw,now,ntw,nnw 2780 2003 format(5i7,2i10) 2781 read(lfntop,2001) cdum 2782 do 105 i=1,naw 2783 read(lfntop,2001) cdum 2784 105 continue 2785 do 106 i=1,nbw*(npars+1) 2786 read(lfntop,2001) cdum 2787 106 continue 2788 do 107 i=1,nhw*(npars+1) 2789 read(lfntop,2001) cdum 2790 107 continue 2791 do 108 i=1,ndw*(npars+1) 2792 read(lfntop,2001) cdum 2793 108 continue 2794 do 109 i=1,now*(npars+1) 2795 read(lfntop,2001) cdum 2796 109 continue 2797 if(ntw.gt.0) then 2798 read(lfntop,2004) (idum,i=1,ntw) 2799 read(lfntop,2004) (idum,i=1,ntw) 2800 2004 format(11i7) 2801 endif 2802 if(nnw.gt.0) then 2803 read(lfntop,2005) (idum,i=1,nnw) 2804 read(lfntop,2005) (idum,i=1,nnw) 2805 2005 format(11i7) 2806 endif 2807 read(lfntop,2001) cdum 2808 do 204 i=1,npars 2809 read(lfntop,2001) cdum 2810 204 continue 2811c 2812 if(nsa.gt.0) then 2813 number=0 2814 ncyc=nsa/nb+1 2815 nums=nb 2816 do 6 icyc=1,ncyc 2817 if(nsa-number.lt.nums) nums=nsa-number 2818 do 7 node=np-1,0,-1 2819 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 2820 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 2821 nsan=ipl(2,2) 2822 if(nsan.gt.0) then 2823 call ga_distribution(ga_is,node,ili,ihi,jli,jhi) 2824 if(npack.eq.0) then 2825 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+lsdyn-1,isl,msa) 2826 else 2827 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa) 2828 call sp_unpack(nsan,isl,islp) 2829 endif 2830 call ga_distribution(ga_s,node,ils,ihs,jls,jhs) 2831 call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa) 2832 call ga_get(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa) 2833 do 8 i=1,nsan 2834 j=isl(i,lsgan)-number 2835 if(j.gt.0.and.j.le.nums) then 2836 bxs(j,1)=xs(i,1) 2837 bxs(j,2)=xs(i,2) 2838 bxs(j,3)=xs(i,3) 2839 bfs(j,1)=fs(i,1) 2840 bfs(j,2)=fs(i,2) 2841 bfs(j,3)=fs(i,3) 2842 ibs(j)=isl(i,lsdyn) 2843 endif 2844 8 continue 2845 endif 2846 7 continue 2847 do 9 i=1,nums 2848 read(lfntop,2009) cat,ism,iss 2849 2009 format(a16,3x,2i7) 2850c 2009 format(a16,3x,2i7) 2851 idyn=iand(ibs(i),3) 2852 if(npener.le.0) then 2853 write(lfnout,1009) ism,iss,cat,idyn, 2854 + (bxs(i,j),j=1,3),(bfs(i,j),j=1,3) 2855 1009 format(2i5,':',a16,i1,1x,3f8.4,3x,3f12.3) 2856 else 2857 write(lfnout,1010) ism,iss,cat,idyn, 2858 + (bxs(i,j),j=1,3),(bfs(i,j),j=1,3),esa(i,1),esa(i,2), 2859 + esa(i,1)+esa(i,2) 2860 1010 format(2i5,':',a16,i1,1x,3f8.4,3x,3f12.3,3(2x,1pe10.3)) 2861 endif 2862 9 continue 2863 6 continue 2864 endif 2865c 2866 close(unit=lfntop,status='keep') 2867c 2868 return 2869 end 2870 subroutine sp_wtrest(lfn) 2871c 2872 implicit none 2873c 2874#include "sp_common.fh" 2875#include "mafdecls.fh" 2876c 2877 integer lfn 2878c 2879 integer i_ltmp,l_ltmp 2880c 2881 if(.not.ma_push_get(mt_int,(mbbl+1)*mbb2,'ltmp',l_ltmp,i_ltmp)) 2882 + call md_abort('Failed to allocate memory for ltmp',0) 2883 call sp_wrest(lfn,int_mb(i_bb),int_mb(i_ltmp),mbbl+1, 2884 + dbl_mb(i_boxs)) 2885 if(.not.ma_pop_stack(l_ltmp)) 2886 + call md_abort('Failed to deallocate ltmp',0) 2887c 2888 return 2889 end 2890 subroutine sp_wrest(lfn,lbbl,ltemp,mdim,boxsiz) 2891c 2892 implicit none 2893c 2894#include "sp_common.fh" 2895#include "msgids.fh" 2896#include "global.fh" 2897c 2898 integer lfn,mdim 2899 integer lbbl(mbbl,mbb2),ltemp(mdim,mbb2) 2900 real*8 boxsiz(maxbox,3) 2901c 2902 integer i,j,k 2903c 2904 if(me.eq.0) then 2905 write(lfn,1000) 2906 1000 format('restart space') 2907 write(lfn,1001) np,mbbl,npx,npy,npz,nbx,nby,nbz,np 2908 1001 format(9i7) 2909 write(lfn,1002) (boxsiz(i,1),i=1,nbx) 2910 write(lfn,1002) (boxsiz(i,2),i=1,nby) 2911 write(lfn,1002) (boxsiz(i,3),i=1,nbz) 2912 1002 format(4e20.12) 2913 endif 2914c 2915 do 1 i=1,np 2916 if(i.eq.me+1) then 2917 do 2 k=1,mbb2 2918 ltemp(1,k)=0 2919 do 3 j=1,nbbl 2920 ltemp(j+1,k)=lbbl(j,k) 2921 3 continue 2922 2 continue 2923 ltemp(1,1)=nbbl 2924 else 2925 do 4 k=1,mbb2 2926 do 5 j=1,mdim 2927 ltemp(j,k)=0 2928 5 continue 2929 4 continue 2930 endif 2931c 2932 call ga_igop(msp_09,ltemp,mdim*mbb2,'+') 2933c 2934 if(me.eq.0) then 2935 write(lfn,1003) i-1,ltemp(1,1) 2936 1003 format(2i7) 2937 do 6 j=1,ltemp(1,1) 2938 write(lfn,1004) (ltemp(j+1,k),k=1,4) 2939 1004 format(8i10) 2940 6 continue 2941 endif 2942c 2943 1 continue 2944c 2945 return 2946 end 2947 logical function sp_rdrest(lfn,fil,boxsiz) 2948c 2949 implicit none 2950c 2951#include "sp_common.fh" 2952#include "mafdecls.fh" 2953c 2954 integer lfn 2955 character*255 fil 2956 real*8 boxsiz(maxbox,3) 2957c 2958 logical sp_rrest 2959 external sp_rrest 2960c 2961 integer i_ltmp,l_ltmp 2962c 2963 if(.not.ma_push_get(mt_int,(mbbl+1)*mbb2,'ltmp',l_ltmp,i_ltmp)) 2964 + call md_abort('Failed to allocate memory for ltmp',0) 2965 sp_rdrest=sp_rrest(lfn,fil,int_mb(i_bb),int_mb(i_ltmp), 2966 + mbbl+1,boxsiz) 2967 if(.not.ma_pop_stack(l_ltmp)) 2968 + call md_abort('Failed to deallocate ltmp',0) 2969c 2970 return 2971 end 2972 logical function sp_rrest(lfn,fil,lbbl,ltemp,mdim,boxsiz) 2973c 2974 implicit none 2975c 2976#include "sp_common.fh" 2977#include "mafdecls.fh" 2978#include "msgids.fh" 2979c 2980 integer lfn,mdim 2981 character*255 fil 2982 integer lbbl(mbbl,mbb2),ltemp(mdim,mbb2) 2983c 2984 integer i,j,k,node,npp,nbl,nbytes,nbxs 2985 integer npxp,npyp,npzp,nbxp,nbyp,nbzp,mbblx 2986 character*13 string 2987 real*8 boxsiz(maxbox,3) 2988c 2989 if(me.eq.0) then 2990 open(unit=lfn,file=fil(1:index(fil,' ')-1), 2991 + status='old',form='formatted',err=9999) 2992 rewind(lfn) 2993c 2994 1 continue 2995 npp=0 2996 read(lfn,1000,end=9997) string 2997 1000 format(a13) 2998 if(string.ne.'restart space') goto 1 2999 read(lfn,1001) npp,mbblx,npxp,npyp,npzp,nbxp,nbyp,nbzp,nbxs 3000 1001 format(9i7) 3001 if(mbblp.lt.mbblx) npp=0 3002 if(npxp.ne.npx) npp=0 3003 if(npyp.ne.npy) npp=0 3004 if(npzp.ne.npz) npp=0 3005 if(nbxp.ne.nbx) npp=0 3006 if(nbyp.ne.nby) npp=0 3007 if(nbzp.ne.nbz) npp=0 3008 if(nbxs.gt.0.and.npp.gt.0) then 3009 read(lfn,1002) (boxsiz(i,1),i=1,nbxp) 3010 read(lfn,1002) (boxsiz(i,2),i=1,nbyp) 3011 read(lfn,1002) (boxsiz(i,3),i=1,nbzp) 3012 1002 format(4e20.12) 3013 endif 3014 9997 continue 3015 do 2 i=1,mbb2 3016 ltemp(1,i)=0 3017 2 continue 3018 endif 3019c 3020 nbytes=ma_sizeof(mt_int,1,mt_byte) 3021 call ga_brdcst(msp_10,npp,nbytes,0) 3022c 3023 if(np.ne.npp) goto 9998 3024c 3025 nbytes=ma_sizeof(mt_dbl,3*maxbox,mt_byte) 3026 call ga_brdcst(msp_10,boxsiz,nbytes,0) 3027c 3028 nbytes=mdim*mbb2*ma_sizeof(mt_int,1,mt_byte) 3029c 3030 do 3 i=1,np 3031c 3032 if(me.eq.0) then 3033 read(lfn,1003) node,nbl 3034 1003 format(2i7) 3035 ltemp(1,1)=node 3036 ltemp(1,2)=nbl 3037 do 4 j=1,nbl 3038 read(lfn,1004) (ltemp(j+1,k),k=1,4) 3039 1004 format(8i10) 3040 if(mbb2.gt.4) then 3041 do 44 k=5,mbb2 3042 ltemp(j+1,k)=0 3043 44 continue 3044 endif 3045 4 continue 3046 endif 3047c 3048 call ga_brdcst(msp_11,ltemp,nbytes,0) 3049c 3050 if(ltemp(1,1).eq.me) then 3051 nbbl=ltemp(1,2) 3052 do 5 k=1,mbb2 3053 do 6 j=1,nbbl 3054 lbbl(j,k)=ltemp(j+1,k) 3055 6 continue 3056 5 continue 3057 endif 3058c 3059 3 continue 3060c 3061 nable=1 3062 sp_rrest=.true. 3063 return 3064 9998 continue 3065 nable=2 3066 sp_rrest=.false. 3067 return 3068 9999 continue 3069 call md_abort('Failed to open restart file',0) 3070 sp_rrest=.false. 3071 return 3072 end 3073 subroutine sp_wrttrj(lfntrj,lxw,lvw,lfw,lxs,lvs,lfs, 3074 + stime,pres,temp,tempw,temps,iwl,xw,vw,fw,xwcr,isl,xs,vs,fs) 3075c 3076 implicit none 3077c 3078#include "sp_common.fh" 3079#include "mafdecls.fh" 3080#include "global.fh" 3081c 3082 integer lfntrj 3083 logical lxw,lvw,lfw,lxs,lvs,lfs 3084 integer iwl(mwm,miw2),isl(msa,mis2) 3085 real*8 stime,pres,temp,tempw,temps 3086 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3) 3087 real*8 xs(msa,3),vs(msa,3),fs(msa,3) 3088c 3089 integer lenscr 3090c 3091 lenscr=ma_inquire_avail(mt_byte)/ 3092 + ((9*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+ 3093 + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1 3094 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx)) 3095 + call md_abort('Failed to allocate bx',0) 3096 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv)) 3097 + call md_abort('Failed to allocate bv',0) 3098 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bf',l_bf,i_bf)) 3099 + call md_abort('Failed to allocate bf',0) 3100 if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br)) 3101 + call md_abort('Failed to allocate br',0) 3102 if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi)) 3103 + call md_abort('Failed to allocate bi',0) 3104 if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n)) 3105 + call md_abort('Failed to allocate n',0) 3106c 3107 call sp_wttrj(lfntrj,lxw,lvw,lfw,lxs,lvs,lfs, 3108 + stime,pres,temp,tempw,temps, 3109 + iwl,int_mb(i_packw),xw,vw,fw,xwcr,isl,int_mb(i_pack),xs,vs,fs, 3110 + int_mb(i_ipl),lenscr,int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv), 3111 + dbl_mb(i_bf),dbl_mb(i_br),int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv), 3112 + dbl_mb(i_bf)) 3113c 3114 if(.not.ma_pop_stack(l_n)) 3115 + call md_abort('Failed to deallocate n',0) 3116 if(.not.ma_pop_stack(l_bi)) 3117 + call md_abort('Failed to deallocate bi',0) 3118 if(.not.ma_pop_stack(l_br)) 3119 + call md_abort('Failed to deallocate br',0) 3120 if(.not.ma_pop_stack(l_bf)) 3121 + call md_abort('Failed to deallocate bf',0) 3122 if(.not.ma_pop_stack(l_bv)) 3123 + call md_abort('Failed to deallocate bv',0) 3124 if(.not.ma_pop_stack(l_bx)) 3125 + call md_abort('Failed to deallocate bx',0) 3126c 3127 return 3128 end 3129 subroutine sp_wttrj(lfntrj,lxw,lvw,lfw,lxs,lvs,lfs, 3130 + stime,pres,temp,tempw,temps,iwl,iwlp,xw,vw,fw,xwcr, 3131 + isl,islp,xs,vs,fs,ipl,nb,ibw,bxw,bvw,bfw,brw,ibs,bxs,bvs,bfs) 3132c 3133 implicit none 3134c 3135#include "sp_common.fh" 3136#include "mafdecls.fh" 3137#include "global.fh" 3138c 3139 integer lfntrj,nb 3140 logical lxw,lvw,lfw,lxs,lvs,lfs 3141 real*8 stime,pres,temp,tempw,temps 3142 integer iwl(mwm,miw2),isl(msa,mis2) 3143 integer iwlp(mwm,npackw),islp(msa,npack) 3144 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3) 3145 real*8 xs(msa,3),vs(msa,3),fs(msa,3) 3146 integer ipl(mbox,mip2),ibw(nb),ibs(nb) 3147 real*8 bxw(nb,3,mwa),bvw(nb,3,mwa),bfw(nb,3,mwa),brw(nb,3) 3148 real*8 bxs(nb,3),bvs(nb,3),bfs(nb,3) 3149c 3150 integer i,j,k,node,ncyc,icyc,numw,nums,number,nwmn,nsan 3151 integer ilp,ihp,jlp,jhp,ili,ihi,jli,jhi,ilw,ihw,jlw,jhw 3152 integer ils,ihs,jls,jhs 3153 character*10 rdate,rtime 3154c 3155 logical lpw,lps 3156c 3157 real*8 dumdst,dist 3158 integer nonh 3159c 3160 dumdst=0.02d0 3161c 3162 lpw=.false. 3163 lps=.false. 3164c 3165 if(me.eq.0) then 3166c 3167 call swatch(rdate,rtime) 3168c 3169 write(lfntrj,1000) 3170 1000 format('frame') 3171 write(lfntrj,1001) stime,temp,pres,rdate,rtime 3172 1001 format(2f12.6,1pe12.5,1x,2a10) 3173 write(lfntrj,1002) ((vlat(i,j),j=1,3),i=1,3) 3174 1002 format(3f12.6) 3175 write(lfntrj,1003) lxw,lvw,lfw,lpw,lxs,lvs,lfs,lps,nwm,nwa,nsa 3176 1003 format(8l1,3i10) 3177c 3178 if((lxw.or.lvw.or.lfw).and.nwm.gt.0) then 3179 number=0 3180 ncyc=nwm/nb+1 3181 numw=nb 3182 do 1 icyc=1,ncyc 3183 if(nwm-number.lt.numw) numw=nwm-number 3184c 3185c begin test code 10/31/2001 3186c initialize ibw to check that all atoms have been received 3187c 3188 do 1112 i=1,nb 3189 ibw(i)=-1 3190 1112 continue 3191c 3192c end test code 3193c 3194 do 2 node=np-1,0,-1 3195 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 3196 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 3197 nwmn=ipl(1,2) 3198 if(nwmn.gt.0) then 3199 call ga_distribution(ga_iw,node,ili,ihi,jli,jhi) 3200 if(npackw.eq.0) then 3201 call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+lwdyn-1,iwl,mwm) 3202 else 3203 call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm) 3204 call sp_unpackw(nwmn,iwl,iwlp) 3205 endif 3206 call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw) 3207 call ga_get(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm) 3208 if(lvw) 3209 + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm) 3210 if(lfw) 3211 + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa+3,jlw+9*mwa+2,fw,mwm) 3212 do 3 i=1,nwmn 3213 j=iwl(i,lwgmn)-number 3214 if(j.gt.0.and.j.le.numw) then 3215 do 4 k=1,nwa 3216 bxw(j,1,k)=xw(i,1,k) 3217 bxw(j,2,k)=xw(i,2,k) 3218 bxw(j,3,k)=xw(i,3,k) 3219 4 continue 3220 if(lvw) then 3221 do 5 k=1,nwa 3222 bvw(j,1,k)=vw(i,1,k) 3223 bvw(j,2,k)=vw(i,2,k) 3224 bvw(j,3,k)=vw(i,3,k) 3225 5 continue 3226 endif 3227 if(lfw) then 3228 do 51 k=1,nwa 3229 bfw(j,1,k)=fw(i,1,k) 3230 bfw(j,2,k)=fw(i,2,k) 3231 bfw(j,3,k)=fw(i,3,k) 3232 51 continue 3233 endif 3234 ibw(j)=iwl(i,lwdyn) 3235 endif 3236 3 continue 3237 endif 3238 2 continue 3239c begin solvent output 3240c low precision 3241 if(nprec.eq.0) then 3242 if(lfw) then 3243 if(lvw) then 3244 do 61 i=1,numw 3245 write(lfntrj,1014) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3), 3246 + (bfw(i,j,k),j=1,3),k=1,nwa) 3247 1014 format(6f8.3,3f8.1) 3248 if(ibw(i).lt.0) call md_abort('Missing solvent in wttrj',i) 3249 61 continue 3250 else 3251 do 71 i=1,numw 3252 write(lfntrj,1015) ((bxw(i,j,k),j=1,3),(bfw(i,j,k),j=1,3),k=1,nwa) 3253 1015 format(3f8.3,3f8.1) 3254 if(ibw(i).lt.0) call md_abort('Missing solvent in wtrst',i) 3255 71 continue 3256 endif 3257 else 3258 if(lvw) then 3259 do 6 i=1,numw 3260 write(lfntrj,1004) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa) 3261 1004 format(6f8.3) 3262 if(ibw(i).lt.0) call md_abort('Missing solvent in wttrj',i) 3263 6 continue 3264 else 3265 do 7 i=1,numw 3266 write(lfntrj,1005) ((bxw(i,j,k),j=1,3),k=1,nwa) 3267 1005 format(3f8.3) 3268 if(ibw(i).lt.0) call md_abort('Missing solvent in wtrst',i) 3269 7 continue 3270 endif 3271 endif 3272c high precision 3273 else 3274 if(lfw) then 3275 if(lvw) then 3276 do 261 i=1,numw 3277 write(lfntrj,2014) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3), 3278 + (bfw(i,j,k),j=1,3),k=1,nwa) 3279 2014 format(6e12.6,/,3e12.6) 3280 if(ibw(i).lt.0) call md_abort('Missing solvent in wttrj',i) 3281 261 continue 3282 else 3283 do 271 i=1,numw 3284 write(lfntrj,2015) ((bxw(i,j,k),j=1,3),(bfw(i,j,k),j=1,3),k=1,nwa) 3285 2015 format(3e12.6,/,3e12.6) 3286 if(ibw(i).lt.0) call md_abort('Missing solvent in wtrst',i) 3287 271 continue 3288 endif 3289 else 3290 if(lvw) then 3291 do 26 i=1,numw 3292 write(lfntrj,2004) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa) 3293 2004 format(6e12.6) 3294 if(ibw(i).lt.0) call md_abort('Missing solvent in wttrj',i) 3295 26 continue 3296 else 3297 do 27 i=1,numw 3298 write(lfntrj,2005) ((bxw(i,j,k),j=1,3),k=1,nwa) 3299 2005 format(3e12.6) 3300 if(ibw(i).lt.0) call md_abort('Missing solvent in wtrst',i) 3301 27 continue 3302 endif 3303 endif 3304 endif 3305c end of solvent output 3306 number=number+numw 3307 1 continue 3308 endif 3309c 3310 if((lxs.or.lvs.or.lfs).and.nsa.gt.0) then 3311 number=0 3312 ncyc=nsa/nb+1 3313 nums=nb 3314 do 8 icyc=1,ncyc 3315 if(nsa-number.lt.nums) nums=nsa-number 3316c 3317c begin test code 10/31/2001 3318c initialize ibw to check that all atoms have been received 3319c 3320 do 1117 i=1,nb 3321 ibs(i)=-1 3322 1117 continue 3323c 3324c end test code 3325c 3326 do 9 node=np-1,0,-1 3327 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 3328 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 3329 nsan=ipl(2,2) 3330 if(nsan.gt.0) then 3331 call ga_distribution(ga_is,node,ili,ihi,jli,jhi) 3332 if(npack.eq.0) then 3333 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+lsdyn-1,isl,msa) 3334 else 3335 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa) 3336 call sp_unpack(nsan,isl,islp) 3337 endif 3338 call ga_distribution(ga_s,node,ils,ihs,jls,jhs) 3339 call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa) 3340 if(lvs) call ga_get(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa) 3341 if(lfs) call ga_get(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa) 3342 nonh=0 3343 do 10 i=1,nsan 3344 j=isl(i,lsgan)-number 3345 if(j.gt.0.and.j.le.nums) then 3346 bxs(j,1)=xs(i,1) 3347 bxs(j,2)=xs(i,2) 3348 bxs(j,3)=xs(i,3) 3349 if(isl(i,lshop).eq.0) nonh=i 3350 if(iand(isl(i,lshop),1).eq.1) then 3351 if(nonh.gt.0) then 3352 bxs(j,1)=xs(nonh,1)-xs(i,1) 3353 bxs(j,2)=xs(nonh,2)-xs(i,2) 3354 bxs(j,3)=xs(nonh,3)-xs(i,3) 3355 dist=sqrt(bxs(j,1)*bxs(j,1)+bxs(j,2)*bxs(j,2)+bxs(j,3)*bxs(j,3)) 3356 dist=dumdst/dist 3357 bxs(j,1)=xs(nonh,1)-dist*bxs(j,1) 3358 bxs(j,2)=xs(nonh,2)-dist*bxs(j,2) 3359 bxs(j,3)=xs(nonh,3)-dist*bxs(j,3) 3360 endif 3361 endif 3362 if(lvs) then 3363 bvs(j,1)=vs(i,1) 3364 bvs(j,2)=vs(i,2) 3365 bvs(j,3)=vs(i,3) 3366 endif 3367 if(lfs) then 3368 bfs(j,1)=fs(i,1) 3369 bfs(j,2)=fs(i,2) 3370 bfs(j,3)=fs(i,3) 3371 endif 3372 ibs(j)=isl(i,lsdyn) 3373 endif 3374 10 continue 3375 endif 3376 9 continue 3377c begin solute output 3378c low precision 3379 if(nprec.eq.0) then 3380 if(lfs) then 3381 if(lvs) then 3382 do 111 i=1,nums 3383 write(lfntrj,1016) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3), 3384 + (bfs(i,j),j=1,3) 3385 1016 format(6f8.3,3f8.1) 3386 if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i) 3387 111 continue 3388 else 3389 do 121 i=1,nums 3390 write(lfntrj,1017) (bxs(i,j),j=1,3),(bfs(i,j),j=1,3) 3391 1017 format(3f8.3,3f8.1) 3392 if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i) 3393 121 continue 3394 endif 3395 else 3396 if(lvs) then 3397 do 11 i=1,nums 3398 write(lfntrj,1006) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3) 3399 1006 format(6f8.3) 3400 if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i) 3401 11 continue 3402 else 3403 do 12 i=1,nums 3404 write(lfntrj,1007) (bxs(i,j),j=1,3) 3405 1007 format(3f8.3) 3406 if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i) 3407 12 continue 3408 endif 3409 endif 3410 else 3411c high precision 3412 if(lfs) then 3413 if(lvs) then 3414 do 2111 i=1,nums 3415 write(lfntrj,2016) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3), 3416 + (bfs(i,j),j=1,3) 3417 2016 format(6e12.6,/,3e12.6) 3418 if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i) 3419 2111 continue 3420 else 3421 do 2121 i=1,nums 3422 write(lfntrj,2017) (bxs(i,j),j=1,3),(bfs(i,j),j=1,3) 3423 2017 format(3e12.6,/,3e12.6) 3424 if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i) 3425 2121 continue 3426 endif 3427 else 3428 if(lvs) then 3429 do 211 i=1,nums 3430 write(lfntrj,2006) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3) 3431 2006 format(6e12.6) 3432 if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i) 3433 211 continue 3434 else 3435 do 212 i=1,nums 3436 write(lfntrj,2007) (bxs(i,j),j=1,3) 3437 2007 format(3e12.6) 3438 if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i) 3439 212 continue 3440 endif 3441 endif 3442 end if 3443c end solute output 3444 number=number+nums 3445 8 continue 3446 endif 3447c 3448 endif 3449c 3450 call util_flush(lfntrj) 3451c 3452 return 3453 end 3454 logical function sp_skip(lfntri,nskip) 3455c 3456 implicit none 3457c 3458#include "sp_common.fh" 3459c 3460 integer lfntri,nskip 3461c 3462 integer i 3463 character*80 card 3464c 3465 if(me.eq.0) then 3466 do 1 i=1,nskip 3467 2 continue 3468 read(lfntri,1000,end=9999) card 3469 1000 format(a) 3470 if(card(1:5).ne.'frame') goto 2 3471 1 continue 3472 endif 3473c 3474 sp_skip=.true. 3475c 3476 return 3477c 3478 9999 continue 3479c 3480 sp_skip=.false. 3481c 3482 return 3483 end 3484 logical function sp_rdtrj(lfntri,lxw,lvw,lfw,lxs,lvs,lfs, 3485 + stime,pres,temp,tempw,temps,iwl,xw,vw,fw,xwcr,isl,xs,vs,fs, 3486 + numwm,numsa) 3487c 3488 implicit none 3489c 3490#include "sp_common.fh" 3491#include "mafdecls.fh" 3492#include "global.fh" 3493c 3494 integer lfntri 3495 logical lxw,lvw,lfw,lxs,lvs,lfs 3496 integer numwm,numsa,iwl(mwm,miw2),isl(msa,mis2) 3497 real*8 stime,pres,temp,tempw,temps 3498 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3) 3499 real*8 xs(msa,3),vs(msa,3),fs(msa,3) 3500c 3501 integer lenscr 3502 logical done 3503c 3504 external sp_rtrj 3505 logical sp_rtrj 3506c 3507 lenscr=ma_inquire_avail(mt_byte)/ 3508 + ((9*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+ 3509 + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1 3510 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx)) 3511 + call md_abort('Failed to allocate bx',0) 3512 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv)) 3513 + call md_abort('Failed to allocate bv',0) 3514 if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bf',l_bf,i_bf)) 3515 + call md_abort('Failed to allocate bf',0) 3516 if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br)) 3517 + call md_abort('Failed to allocate br',0) 3518 if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi)) 3519 + call md_abort('Failed to allocate bi',0) 3520 if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n)) 3521 + call md_abort('Failed to allocate n',0) 3522c 3523 done=sp_rtrj(lfntri,lxw,lvw,lfw,lxs,lvs,lfs, 3524 + stime,pres,temp,tempw,temps, 3525 + iwl,int_mb(i_packw),xw,vw,fw,xwcr,isl,int_mb(i_pack),xs,vs,fs, 3526 + int_mb(i_ipl),lenscr,int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv), 3527 + dbl_mb(i_bf),dbl_mb(i_br),int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv), 3528 + dbl_mb(i_bf)) 3529c 3530 if(.not.ma_pop_stack(l_n)) 3531 + call md_abort('Failed to deallocate n',0) 3532 if(.not.ma_pop_stack(l_bi)) 3533 + call md_abort('Failed to deallocate bi',0) 3534 if(.not.ma_pop_stack(l_br)) 3535 + call md_abort('Failed to deallocate br',0) 3536 if(.not.ma_pop_stack(l_bf)) 3537 + call md_abort('Failed to deallocate bf',0) 3538 if(.not.ma_pop_stack(l_bv)) 3539 + call md_abort('Failed to deallocate bv',0) 3540 if(.not.ma_pop_stack(l_bx)) 3541 + call md_abort('Failed to deallocate bx',0) 3542c 3543 call ga_sync() 3544c 3545 call sp_gagetixv(me,iwl,int_mb(i_packw),xw,xwcr,vw,numwm, 3546 + isl,int_mb(i_pack),xs,vs,numsa,int_mb(i_ipl)) 3547c 3548 sp_rdtrj=done 3549c 3550 return 3551 end 3552 subroutine sp_gettp(temp,pres) 3553c 3554 implicit none 3555c 3556#include "sp_common.fh" 3557c 3558 real*8 temp,pres 3559c 3560 temp=tmp 3561 pres=prs 3562c 3563 return 3564 end 3565c 3566 logical function sp_rtrj(lfntri,lxw,lvw,lfw,lxs,lvs,lfs, 3567 + stime,pres,temp,tempw,temps,iwl,iwlp,xw,vw,fw,xwcr, 3568 + isl,islp,xs,vs,fs,ipl,nb,ibw,bxw,bvw,bfw,brw,ibs,bxs,bvs,bfs) 3569c 3570 implicit none 3571c 3572#include "sp_common.fh" 3573#include "mafdecls.fh" 3574#include "global.fh" 3575#include "msgids.fh" 3576c 3577 integer lfntri,nb 3578 logical lxw,lvw,lfw,lxs,lvs,lfs 3579 real*8 stime,pres,temp,tempw,temps 3580 integer iwl(mwm,miw2),isl(msa,mis2) 3581 integer iwlp(mwm,npackw),islp(msa,npack) 3582 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3) 3583 real*8 xs(msa,3),vs(msa,3),fs(msa,3) 3584 integer ipl(mbox,mip2),ibw(nb),ibs(nb) 3585 real*8 bxw(nb,3,mwa),bvw(nb,3,mwa),bfw(nb,3,mwa),brw(nb,3) 3586 real*8 bxs(nb,3),bvs(nb,3),bfs(nb,3) 3587c 3588 integer i,j,k,node,ncyc,icyc,numw,nums,number,nwmn,nsan 3589 integer ilp,ihp,jlp,jhp,ili,ihi,jli,jhi,ilw,ihw,jlw,jhw 3590 integer ils,ihs,jls,jhs 3591 character*10 rdate,rtime 3592 character*80 card 3593c 3594 logical lpw,lps 3595c 3596 lpw=.false. 3597 lps=.false. 3598c 3599 if(me.eq.0) then 3600c 3601 100 continue 3602 read(lfntri,1000,end=9999) card 3603 1000 format(a) 3604 if(card(1:5).ne.'frame') goto 100 3605c 3606 read(lfntri,1001) stime,temp,pres,rdate,rtime 3607 1001 format(2f12.6,1pe12.5,1x,2a10) 3608 tmp=temp 3609 prs=pres 3610 read(lfntri,1002) ((vlat(i,j),j=1,3),i=1,3) 3611 1002 format(3f12.6) 3612 read(lfntri,1000,end=9999) card 3613c 3614 if(card(8:8).eq.'T'.or.card(8:8).eq.'F') then 3615 read(card,1003) lxw,lvw,lfw,lpw,lxs,lvs,lfs,lps,nwm,nwa,nsa 3616 1003 format(8l1,3i10) 3617 elseif(card(6:6).eq.'T'.or.card(6:6).eq.'F') then 3618 read(card,1023) lxw,lvw,lfw,lxs,lvs,lfs,nwm,nwa,nsa 3619 1023 format(6l1,3i10) 3620 lpw=.false. 3621 lps=.false. 3622 else 3623 read(card,1033) lxw,lvw,lxs,lvs,nwm,nwa,nsa 3624 1033 format(4l1,3i10) 3625 lpw=.false. 3626 lps=.false. 3627 lfw=.false. 3628 lfs=.false. 3629 endif 3630c 3631 if((lxw.or.lvw.or.lfw).and.nwm.gt.0) then 3632 number=0 3633 ncyc=nwm/nb+1 3634 numw=nb 3635 do 1 icyc=1,ncyc 3636 if(nwm-number.lt.numw) numw=nwm-number 3637c 3638 if(nprec.eq.0) then 3639 if(lfw) then 3640 if(lvw) then 3641 do 61 i=1,numw 3642 read(lfntri,1014) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3), 3643 + (bfw(i,j,k),j=1,3),k=1,nwa) 3644 1014 format(6f8.3,3f8.1) 3645 61 continue 3646 else 3647 do 71 i=1,numw 3648 read(lfntri,1015) ((bxw(i,j,k),j=1,3),(bfw(i,j,k),j=1,3),k=1,nwa) 3649 1015 format(3f8.3,3f8.1) 3650 71 continue 3651 endif 3652 else 3653 if(lvw) then 3654 do 6 i=1,numw 3655 read(lfntri,1004) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa) 3656 1004 format(6f8.3) 3657 6 continue 3658 else 3659 do 7 i=1,numw 3660 read(lfntri,1005) ((bxw(i,j,k),j=1,3),k=1,nwa) 3661 1005 format(3f8.3) 3662 7 continue 3663 endif 3664 endif 3665 else 3666 if(lfw) then 3667 if(lvw) then 3668 do 261 i=1,numw 3669 read(lfntri,2014) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3), 3670 + (bfw(i,j,k),j=1,3),k=1,nwa) 3671 2014 format(6e12.6,/,3e12.6) 3672 261 continue 3673 else 3674 do 271 i=1,numw 3675 read(lfntri,2015) ((bxw(i,j,k),j=1,3),(bfw(i,j,k),j=1,3),k=1,nwa) 3676 2015 format(3e12.6,/,3e12.6) 3677 271 continue 3678 endif 3679 else 3680 if(lvw) then 3681 do 26 i=1,numw 3682 read(lfntri,2004) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa) 3683 2004 format(6e12.6) 3684 26 continue 3685 else 3686 do 27 i=1,numw 3687 read(lfntri,2005) ((bxw(i,j,k),j=1,3),k=1,nwa) 3688 2005 format(3e12.6) 3689 27 continue 3690 endif 3691 endif 3692 endif 3693c 3694 do 2 node=np-1,0,-1 3695 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 3696 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 3697 nwmn=ipl(1,2) 3698 if(nwmn.gt.0) then 3699 call ga_distribution(ga_iw,node,ili,ihi,jli,jhi) 3700 if(npackw.eq.0) then 3701 call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+lwdyn-1,iwl,mwm) 3702 else 3703 call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm) 3704 call sp_unpackw(nwmn,iwl,iwlp) 3705 endif 3706 call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw) 3707 call ga_get(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm) 3708 if(lvw) 3709 + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm) 3710 if(lfw) 3711 + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa+3,jlw+9*mwa+2,fw,mwm) 3712 do 3 i=1,nwmn 3713 j=iwl(i,lwgmn)-number 3714 if(j.gt.0.and.j.le.numw) then 3715 do 4 k=1,nwa 3716 xw(i,1,k)=bxw(j,1,k) 3717 xw(i,2,k)=bxw(j,2,k) 3718 xw(i,3,k)=bxw(j,3,k) 3719 4 continue 3720 if(lvw) then 3721 do 5 k=1,nwa 3722 vw(i,1,k)=bvw(j,1,k) 3723 vw(i,2,k)=bvw(j,2,k) 3724 vw(i,3,k)=bvw(j,3,k) 3725 5 continue 3726 endif 3727 if(lfw) then 3728 do 51 k=1,nwa 3729 fw(i,1,k)=bfw(j,1,k) 3730 fw(i,2,k)=bfw(j,2,k) 3731 fw(i,3,k)=bfw(j,3,k) 3732 51 continue 3733 endif 3734 endif 3735 3 continue 3736 call ga_put(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm) 3737 if(lvw) 3738 + call ga_put(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm) 3739 if(lfw) 3740 + call ga_put(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa+3,jlw+9*mwa+2,fw,mwm) 3741 endif 3742 2 continue 3743 number=number+numw 3744 1 continue 3745 endif 3746c 3747 if((lxs.or.lvs.or.lfs).and.nsa.gt.0) then 3748 number=0 3749 ncyc=nsa/nb+1 3750 nums=nb 3751 do 8 icyc=1,ncyc 3752 if(nsa-number.lt.nums) nums=nsa-number 3753c 3754 if(nprec.eq.0) then 3755 if(lfs) then 3756 if(lvs) then 3757 do 111 i=1,nums 3758 read(lfntri,1016) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3), 3759 + (bfs(i,j),j=1,3) 3760 1016 format(6f8.3,3f8.1) 3761 111 continue 3762 else 3763 do 121 i=1,nums 3764 read(lfntri,1017) (bxs(i,j),j=1,3),(bfs(i,j),j=1,3) 3765 1017 format(3f8.3,3f8.1) 3766 121 continue 3767 endif 3768 else 3769 if(lvs) then 3770 do 11 i=1,nums 3771 read(lfntri,1006) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3) 3772 1006 format(6f8.3) 3773 11 continue 3774 else 3775 do 12 i=1,nums 3776 read(lfntri,1007) (bxs(i,j),j=1,3) 3777 1007 format(3f8.3) 3778 12 continue 3779 endif 3780 endif 3781 else 3782 if(lfs) then 3783 if(lvs) then 3784 do 2111 i=1,nums 3785 read(lfntri,2016) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3), 3786 + (bfs(i,j),j=1,3) 3787 2016 format(6e12.6,/,3e12.6) 3788 2111 continue 3789 else 3790 do 2121 i=1,nums 3791 read(lfntri,2017) (bxs(i,j),j=1,3),(bfs(i,j),j=1,3) 3792 2017 format(3e12.6,/,3e12.6) 3793 2121 continue 3794 endif 3795 else 3796 if(lvs) then 3797 do 211 i=1,nums 3798 read(lfntri,2006) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3) 3799 2006 format(6e12.6) 3800 211 continue 3801 else 3802 do 212 i=1,nums 3803 read(lfntri,2007) (bxs(i,j),j=1,3) 3804 2007 format(3e12.6) 3805 212 continue 3806 endif 3807 endif 3808 endif 3809c 3810 do 9 node=np-1,0,-1 3811 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 3812 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 3813 nsan=ipl(2,2) 3814 if(nsan.gt.0) then 3815 call ga_distribution(ga_is,node,ili,ihi,jli,jhi) 3816 if(npack.eq.0) then 3817 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+lsdyn-1,isl,msa) 3818 else 3819 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa) 3820 call sp_unpack(nsan,isl,islp) 3821 endif 3822 call ga_distribution(ga_s,node,ils,ihs,jls,jhs) 3823 call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa) 3824 if(lvs) call ga_get(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa) 3825 if(lfs) call ga_get(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa) 3826 do 10 i=1,nsan 3827 j=isl(i,lsgan)-number 3828 if(j.gt.0.and.j.le.nums) then 3829 xs(i,1)=bxs(j,1) 3830 xs(i,2)=bxs(j,2) 3831 xs(i,3)=bxs(j,3) 3832 if(lvs) then 3833 vs(i,1)=bvs(j,1) 3834 vs(i,2)=bvs(j,2) 3835 vs(i,3)=bvs(j,3) 3836 endif 3837 if(lfs) then 3838 fs(i,1)=bfs(j,1) 3839 fs(i,2)=bfs(j,2) 3840 fs(i,3)=bfs(j,3) 3841 endif 3842 endif 3843 10 continue 3844 call ga_put(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa) 3845 if(lvs) call ga_put(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa) 3846 if(lfs) call ga_put(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa) 3847 endif 3848 9 continue 3849 number=number+nums 3850 8 continue 3851 endif 3852c 3853 endif 3854c 3855 call ga_brdcst(msp_15,vlat, 3856 + 9*ma_sizeof(mt_dbl,1,mt_byte),0) 3857c 3858 box(1)=vlat(1,1) 3859 box(2)=vlat(2,2) 3860 box(3)=vlat(3,3) 3861 boxh(1)=5.0d-1*box(1) 3862 boxh(2)=5.0d-1*box(2) 3863 boxh(3)=5.0d-1*box(3) 3864c 3865 sp_rtrj=.true. 3866c 3867 return 3868c 3869 9999 continue 3870c 3871 sp_rtrj=.false. 3872c 3873 return 3874c 3875 end 3876c 3877 subroutine sp_wrtmro(lfnmro,stime,pres,temp,tempw,temps, 3878 + iwl,xw,vw,xwcr,isl,xs,vs,prjct) 3879c 3880 implicit none 3881c 3882#include "sp_common.fh" 3883#include "mafdecls.fh" 3884#include "global.fh" 3885c 3886 integer lfnmro 3887 integer iwl(mwm,miw2),isl(msa,mis2) 3888 real*8 stime,pres,temp,tempw,temps 3889 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3) 3890 real*8 xs(msa,3),vs(msa,3) 3891 character*80 prjct 3892c 3893 project=prjct 3894c 3895 call sp_wtmro(lfnmro,stime,pres,temp,tempw,temps, 3896 + iwl,int_mb(i_packw),xw,vw,xwcr,isl,int_mb(i_pack),xs,vs, 3897 + int_mb(i_ipl)) 3898c 3899 return 3900 end 3901 subroutine sp_wtmro(lfnmro,stime,pres,temp,tempw,temps, 3902 + iwl,iwlp,xw,vw,xwcr,isl,islp,xs,vs,ipl) 3903c 3904 implicit none 3905c 3906#include "sp_common.fh" 3907#include "mafdecls.fh" 3908#include "global.fh" 3909c 3910 integer lfnmro 3911 real*8 stime,pres,temp,tempw,temps 3912 integer iwl(mwm,miw2),isl(msa,mis2) 3913 integer iwlp(mwm,npackw),islp(msa,npack) 3914 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3) 3915 real*8 xs(msa,3),vs(msa,3) 3916 integer ipl(mbox,mip2) 3917c 3918 integer j,k,l,nwmn,nsan,node,ilp,ihp,jlp,jhp,ili,ihi,jli,jhi 3919 integer ilw,ihw,jlw,jhw,ils,ihs,jls,jhs 3920 character*10 rdate,rtime 3921 character*18 user 3922#ifdef USE_POSIXF 3923 integer ilen,ierror 3924#endif 3925c 3926 write(lfnmro) nwm,nwa,nsa,stime,temp,pres,vlat,nhist 3927c 3928 call swatch(rdate,rtime) 3929#ifdef USE_POSIXF 3930 call pxfgetlogin(user, ilen, ierror) 3931#else 3932 call getlog(user) 3933#endif 3934 if(user(18:18).ne.' ') user=' ' 3935 hist(nhist)(1:18)=user 3936 hist(nhist)(19:28)=rdate 3937 hist(nhist)(29:48)=rtime 3938 hist(nhist)(49:108)=project(1:60) 3939 write(lfnmro) (hist(j),j=1,nhist) 3940c 3941 do 1 node=np-1,0,-1 3942 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 3943 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 3944 write(lfnmro) ((ipl(j,k),j=1,mbox),k=1,mip2) 3945 nwmn=ipl(1,2) 3946 nsan=ipl(2,2) 3947 if(nwmn.gt.0) then 3948 call ga_distribution(ga_iw,node,ili,ihi,jli,jhi) 3949 if(npackw.eq.0) then 3950 call ga_get(ga_iw,ili,ili+nwmn-1,jli,jhi,iwl,mwm) 3951 else 3952 call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm) 3953 call sp_unpackw(nwmn,iwl,iwlp) 3954 endif 3955 call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw) 3956 call ga_get(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm) 3957 call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm) 3958 write(lfnmro) ((iwl(j,k),j=1,nwmn),k=1,miw2) 3959 write(lfnmro) (((xw(j,k,l),j=1,nwmn),k=1,3),l=1,nwa) 3960 write(lfnmro) (((vw(j,k,l),j=1,nwmn),k=1,3),l=1,nwa) 3961 endif 3962 if(nsan.gt.0) then 3963 call ga_distribution(ga_is,node,ili,ihi,jli,jhi) 3964 if(npack.eq.0) then 3965 call ga_get(ga_is,ili,ili+nsan-1,jli,jhi,isl,msa) 3966 else 3967 call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa) 3968 call sp_unpack(nsan,isl,islp) 3969 endif 3970 call ga_distribution(ga_s,node,ils,ihs,jls,jhs) 3971 call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa) 3972 call ga_get(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa) 3973 write(lfnmro) ((isl(j,k),j=1,nsan),k=1,mis2) 3974 write(lfnmro) ((xs(j,k),j=1,nsan),k=1,3) 3975 write(lfnmro) ((vs(j,k),j=1,nsan),k=1,3) 3976 endif 3977 1 continue 3978c 3979 return 3980 end 3981 logical function sp_rdmri(lfnmri,stime,pres,temp,tempw,temps, 3982 + iwl,xw,vw,xwcr,isl,xs,vs) 3983c 3984 implicit none 3985c 3986#include "sp_common.fh" 3987#include "mafdecls.fh" 3988#include "global.fh" 3989c 3990 logical sp_rmri 3991 external sp_rmri 3992c 3993 integer lfnmri 3994 integer iwl(mwm,miw2),isl(msa,mis2) 3995 real*8 stime,pres,temp,tempw,temps 3996 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3) 3997 real*8 xs(msa,3),vs(msa,3) 3998c 3999 sp_rdmri=sp_rmri(lfnmri,stime,pres,temp,tempw,temps, 4000 + iwl,int_mb(i_packw),xw,vw,xwcr,isl,int_mb(i_pack),xs,vs, 4001 + int_mb(i_ipl)) 4002c 4003 return 4004 end 4005 logical function sp_rmri(lfnmri,stime,pres,temp,tempw,temps, 4006 + iwl,iwlp,xw,vw,xwcr,isl,islp,xs,vs,ipl) 4007c 4008 implicit none 4009c 4010#include "sp_common.fh" 4011#include "mafdecls.fh" 4012#include "msgids.fh" 4013#include "global.fh" 4014c 4015 integer lfnmri 4016 real*8 stime,pres,temp,tempw,temps 4017 integer iwl(mwm,miw2),isl(msa,mis2) 4018 integer iwlp(mwm,npackw),islp(msa,npack) 4019 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3) 4020 real*8 xs(msa,3),vs(msa,3) 4021 integer ipl(mbox,mip2) 4022c 4023 integer i,j,k,l,nwmn,nsan,node,ilp,ihp,jlp,jhp,ili,ihi,jli,jhi 4024 integer ilw,ihw,jlw,jhw,ils,ihs,jls,jhs 4025 integer ltemp(3) 4026 real*8 rtemp(12) 4027c 4028 if(me.eq.0) then 4029 read(lfnmri,err=9,end=9) ltemp,rtemp,nhist 4030 if(nhist.gt.0) read(lfnmri,err=9,end=9) (hist(j),j=1,nhist) 4031 do 1 node=np-1,0,-1 4032 read(lfnmri,err=9) ((ipl(j,k),j=1,mbox),k=1,mip2) 4033 nwmn=ipl(1,2) 4034 nsan=ipl(2,2) 4035 call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp) 4036 call ga_put(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 4037 if(nwmn.gt.0) then 4038 read(lfnmri,err=9) ((iwl(j,k),j=1,nwmn),k=1,miw2) 4039 read(lfnmri,err=9) (((xw(j,k,l),j=1,nwmn),k=1,3),l=1,nwa) 4040 read(lfnmri,err=9) (((vw(j,k,l),j=1,nwmn),k=1,3),l=1,nwa) 4041 call ga_distribution(ga_iw,node,ili,ihi,jli,jhi) 4042 if(npackw.eq.0) then 4043 call ga_put(ga_iw,ili,ili+nwmn-1,jli,jhi,iwl,mwm) 4044 else 4045 call sp_packw(nwmn,iwl,iwlp) 4046 call ga_put(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm) 4047 endif 4048 call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw) 4049 call ga_put(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm) 4050 call ga_put(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm) 4051 endif 4052 if(nsan.gt.0) then 4053 read(lfnmri,err=9) ((isl(j,k),j=1,nsan),k=1,mis2) 4054 read(lfnmri,err=9) ((xs(j,k),j=1,nsan),k=1,3) 4055 read(lfnmri,err=9) ((vs(j,k),j=1,nsan),k=1,3) 4056 call ga_distribution(ga_is,node,ili,ihi,jli,jhi) 4057 if(npack.eq.0) then 4058 call ga_put(ga_is,ili,ili+nsan-1,jli,jhi,isl,msa) 4059 else 4060 call sp_pack(nsan,isl,islp) 4061 call ga_put(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa) 4062 endif 4063 call ga_distribution(ga_s,node,ils,ihs,jls,jhs) 4064 call ga_put(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa) 4065 call ga_put(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa) 4066 endif 4067 1 continue 4068 endif 4069c 4070 call ga_brdcst(msp_12,ltemp,3*ma_sizeof(mt_int,1,mt_byte),0) 4071 call ga_brdcst(msp_13,rtemp,12*ma_sizeof(mt_dbl,1,mt_byte),0) 4072 nwm=ltemp(1) 4073 nwa=ltemp(2) 4074 nsa=ltemp(3) 4075 stime=rtemp(1) 4076 temp=rtemp(2) 4077 pres=rtemp(3) 4078 k=3 4079 do 2 i=1,3 4080 do 3 j=1,3 4081 k=k+1 4082 vlat(i,j)=rtemp(k) 4083 3 continue 4084 2 continue 4085 call ga_distribution(ga_ip,me,ilp,ihp,jlp,jhp) 4086 call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox) 4087 nwmloc=ipl(1,2) 4088 nsaloc=ipl(2,2) 4089 if(nwmloc.gt.0) then 4090 call ga_distribution(ga_iw,me,ili,ihi,jli,jhi) 4091 if(npackw.eq.0) then 4092 call ga_get(ga_iw,ili,ili+nwmloc-1,jli,jhi,iwl,mwm) 4093 else 4094 call ga_get(ga_iw,ili,ili+nwmloc-1,jli,jli+npackw-1,iwlp,mwm) 4095 call sp_unpackw(nwmloc,iwl,iwlp) 4096 endif 4097 call ga_distribution(ga_w,me,ilw,ihw,jlw,jhw) 4098 call ga_get(ga_w,ilw,ilw+nwmloc-1,jlw,jlw+3*mwa-1,xw,mwm) 4099 call ga_get(ga_w,ilw,ilw+nwmloc-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm) 4100 endif 4101 if(nsaloc.gt.0) then 4102 call ga_distribution(ga_is,me,ili,ihi,jli,jhi) 4103 if(npack.eq.0) then 4104 call ga_get(ga_is,ili,ili+nsaloc-1,jli,jhi,isl,msa) 4105 else 4106 call ga_get(ga_is,ili,ili+nsaloc-1,jli,jli+npack-1,islp,msa) 4107 call sp_unpack(nsaloc,isl,islp) 4108 endif 4109 call ga_distribution(ga_s,me,ils,ihs,jls,jhs) 4110 call ga_get(ga_s,ils,ils+nsaloc-1,jls,jls+2,xs,msa) 4111 call ga_get(ga_s,ils,ils+nsaloc-1,jls+3,jls+5,vs,msa) 4112 endif 4113c 4114 sp_rmri=.true. 4115 return 4116c 4117 9 continue 4118 sp_rmri=.false. 4119 return 4120 end 4121 subroutine sp_fix() 4122c 4123 implicit none 4124c 4125 return 4126 end 4127