text:
      SUBROUTINE woutp(I,NHW,NHE,NVN,NVS)
C
C     THIS SUBROUTINE WRITES OUT ALL PLANT VARIABLES TO
C     CHECKPOINT FILES AT THE FREQUENCY GIVEN IN THE OPTIONS
C     FILE SO THAT OUTPUTS FROM EARLIER MODEL RUNS CAN BE USED
C     TO INITIALIZE LATER MODEL RUNS
C
      include "parameters.h"
      include "files.h"
      include "blkc.h"
      include "blk1cp.h"
      include "blk1cr.h"
      include "blk1g.h"
      include "blk1n.h"
      include "blk1p.h"
      include "blk1s.h"
      include "blk3.h"
      include "blk5.h"
      include "blk8a.h"
      include "blk8b.h"
      include "blk9a.h"
      include "blk9b.h"
      include "blk9c.h"
      include "blk12a.h"
      include "blk12b.h"
      include "blk14.h"
      include "blk18a.h"
      include "blk18b.h"
      DO 9990 NX=NHW,NHE
      DO 9985 NY=NVN,NVS
      WRITE(26,90)I,IDATA(3),NP(NY,NX),IFLGT(NY,NX)
     2,(IFLGC(NZ,NY,NX),NZ=1,NP(NY,NX))
     3,(IDTH(NZ,NY,NX),NZ=1,NP(NY,NX))
      WRITE(26,97)I,IDATA(3),NP(NY,NX)
     2,((WTSTDG(M,NZ,NY,NX),M=1,4),NZ=1,NP(NY,NX))
     3,((WTSTDN(M,NZ,NY,NX),M=1,4),NZ=1,NP(NY,NX))
     4,((WTSTDP(M,NZ,NY,NX),M=1,4),NZ=1,NP(NY,NX))
      IF(IFLGT(NY,NX).GT.0)THEN
      DO 9980 NZ=1,NP(NY,NX)
      IF(IFLGC(NZ,NY,NX).NE.0)THEN
      WRITE(26,91)I,IDATA(3),NZ,IYR0(NZ,NY,NX),IDAY0(NZ,NY,NX)
     2,IYRH(NZ,NY,NX),IDAYH(NZ,NY,NX),NG(NZ,NY,NX),IDTHP(NZ,NY,NX)
     3,IDTHR(NZ,NY,NX),NBR(NZ,NY,NX),NBT(NZ,NY,NX),NB1(NZ,NY,NX)
     4,IFLGI(NZ,NY,NX),NRT(NZ,NY,NX),NIX(NZ,NY,NX),MY(NZ,NY,NX)
     5,(NINR(NR,NZ,NY,NX),NR=1,NRT(NZ,NY,NX))
      WRITE(26,93)I,IDATA(3),NZ,TCUPTK(NZ,NY,NX),TCSNC(NZ,NY,NX)
     2,TZUPTK(NZ,NY,NX),TZSNC(NZ,NY,NX),TPUPTK(NZ,NY,NX),TPSNC(NZ,NY,NX)
     3,TZUPFX(NZ,NY,NX),TCO2T(NZ,NY,NX),CTRAN(NZ,NY,NX),CARBN(NZ,NY,NX)
     4,TCC(NZ,NY,NX),TKC(NZ,NY,NX),TCG(NZ,NY,NX),TKG(NZ,NY,NX)
     5,TFN3(NZ,NY,NX),WVSTK(NZ,NY,NX),VOLWP(NZ,NY,NX)
     6,PSILT(NZ,NY,NX),PSILO(NZ,NY,NX),PSILG(NZ,NY,NX),WTRTA(NZ,NY,NX)
     7,ARLFP(NZ,NY,NX),ARSTP(NZ,NY,NX),PP(NZ,NY,NX),HVSTC(NZ,NY,NX)
     8,HVSTN(NZ,NY,NX),HVSTP(NZ,NY,NX),THVSTC(NZ,NY,NX),THVSTN(NZ,NY,NX)
     3,THVSTP(NZ,NY,NX),TCO2A(NZ,NY,NX),RSETC(NZ,NY,NX),RSETN(NZ,NY,NX)
     9,RSETP(NZ,NY,NX),TNH3C(NZ,NY,NX),TCSN0(NZ,NY,NX),PPI(NZ,NY,NX)
     1,PPX(NZ,NY,NX),WTSTG(NZ,NY,NX),WTSTGN(NZ,NY,NX),WTSTGP(NZ,NY,NX)
     2,ZC(NZ,NY,NX),VOLWC(NZ,NY,NX),CF(NZ,NY,NX),VCO2F(NZ,NY,NX)
     3,VCH4F(NZ,NY,NX),VOXYF(NZ,NY,NX),VNH3F(NZ,NY,NX),VN2OF(NZ,NY,NX)
     4,VPO4F(NZ,NY,NX)
      WRITE(26,92)I,IDATA(3),NZ,(GROUP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(PSTG(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(PSTGI(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(PSTGF(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(VSTG(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(VSTGX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(GSTGI(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(GSTGF(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(TGSTGI(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(TGSTGF(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(VRNS(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(VRNF(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(VRNY(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(VRNZ(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(ATRP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(FLG4(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,92)I,IDATA(3),NZ,(FLGZ(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(IFLGA(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(IFLGE(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(IFLGF(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(IFLGR(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(IFLGQ(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(IFLGD(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(IDTHB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(NBTB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(KLEAF(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(KVSTG(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(26,95)I,IDATA(3),NZ,(KVSTGN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      DO 9955 N=1,10
9955  WRITE(26,95)I,IDATA(3),NZ,(IDAY(N,NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,96)I,IDATA(3),NZ,CPOOLP(NZ,NY,NX),CPOLNP(NZ,NY,NX)
     3,VHCPC(NZ,NY,NX),DTKC(NZ,NY,NX),WTSHT(NZ,NY,NX)
     4,WTLS(NZ,NY,NX),WTRT(NZ,NY,NX),WTSHN(NZ,NY,NX),WTSHP(NZ,NY,NX)
     5,WTLF(NZ,NY,NX),WTSHE(NZ,NY,NX),WTSTK(NZ,NY,NX),WTRSV(NZ,NY,NX)
     6,WTHSK(NZ,NY,NX),WTEAR(NZ,NY,NX),WTGR(NZ,NY,NX),WTND(NZ,NY,NX)
     7,WTRVC(NZ,NY,NX),WTRVN(NZ,NY,NX),WTRVP(NZ,NY,NX),HTCTL(NZ,NY,NX)
     8,SDPTH(NZ,NY,NX),WSTR(NZ,NY,NX)
     9,CHILL(NZ,NY,NX),WTRTS(NZ,NY,NX),FRADP(NZ,NY,NX)
      WRITE(27,92)I,IDATA(3),NZ,(CPOOL(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(ZPOOL(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(PPOOL(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(CPOLNB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(ZPOLNB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(PPOLNB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSHTB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTLFB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTNDB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSHEB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSTKB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WVSTKB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTRSVB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTHSKB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTEARB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTGRB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTLSB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSHTN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTLFBN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTNDBN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSHBN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSTBN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTRSBN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTHSBN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTEABN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTGRBN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSHTP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTLFBP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTNDBP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSHBP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSTBP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTRSBP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTHSBP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTEABP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTGRBP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(GRNXB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(GRNOB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(GRWTB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(ARLFB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WGLFX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WGLFNX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WGLFPX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(ARLFZ(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(RCZLX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(RCPLX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(RCCLX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WGSHEX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WGSHNX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WGSHPX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(HTSHEX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(RCZSX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(RCPSX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(RCCSX(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSTXB(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSTXN(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,92)I,IDATA(3),NZ,(WTSTXP(NB,NZ,NY,NX),NB=1,NBR(NZ,NY,NX))
      WRITE(27,94)I,IDATA(3),NZ,(ARLFV(L,NZ,NY,NX),L=1,JC)
      WRITE(27,94)I,IDATA(3),NZ,(ARSTV(L,NZ,NY,NX),L=1,JC)
      DO 9945 NB=1,NBR(NZ,NY,NX)
      WRITE(28,93)I,IDATA(3),NZ,(CPOOL4(K,NB,NZ,NY,NX),K=1,25)
      WRITE(28,93)I,IDATA(3),NZ,(CPOOL3(K,NB,NZ,NY,NX),K=1,25)
      WRITE(28,93)I,IDATA(3),NZ,(CO2B(K,NB,NZ,NY,NX),K=1,25)
      WRITE(28,93)I,IDATA(3),NZ,(HCOB(K,NB,NZ,NY,NX),K=1,25)
      WRITE(28,93)I,IDATA(3),NZ,(ARLF(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGLF(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WSLF(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(HTSHE(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGSHE(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WSSHE(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(HTNODE(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(HTNODX(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGNODE(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGLFN(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGSHN(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGNODN(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGLFP(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGSHP(K,NB,NZ,NY,NX),K=0,25)
      WRITE(28,93)I,IDATA(3),NZ,(WGNODP(K,NB,NZ,NY,NX),K=0,25)
      DO 9950 K=0,25
      WRITE(28,94)I,IDATA(3),NZ,(ARLFL(L,K,NB,NZ,NY,NX),L=1,JC)
      WRITE(28,94)I,IDATA(3),NZ,(WGLFL(L,K,NB,NZ,NY,NX),L=1,JC)
      WRITE(28,94)I,IDATA(3),NZ,(WGLFLN(L,K,NB,NZ,NY,NX),L=1,JC)
      WRITE(28,94)I,IDATA(3),NZ,(WGLFLP(L,K,NB,NZ,NY,NX),L=1,JC)
      IF(K.NE.0)THEN
      DO 9940 N=1,4
      WRITE(28,94)I,IDATA(3),NZ,(SURF(N,L,K,NB,NZ,NY,NX),L=1,JC)
9940  CONTINUE
      ENDIF
9950  CONTINUE
      WRITE(28,94)I,IDATA(3),NZ,(ARSTK(L,NB,NZ,NY,NX),L=1,JC)
      DO 9945 L=1,JC
      WRITE(28,92)I,IDATA(3),NZ,(SURFB(N,L,NB,NZ,NY,NX),N=1,4)
9945  CONTINUE
      DO 9970 N=1,MY(NZ,NY,NX)
      WRITE(29,94)I,IDATA(3),NZ,(PSIRT(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(PSIRO(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(PSIRG(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTN1(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTNL(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTLGP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTDNP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTVLP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTVLW(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RRAD1(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RRAD2(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTARP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTLGA(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RCO2M(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RCO2A(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RCO2N(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(CO2A(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(OXYA(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(CH4A(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(Z2OA(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(ZH3A(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(CO2P(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(OXYP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(CH4P(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(Z2OP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(ZH3P(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTRTL(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTRTD(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WSRTL(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(ROXYP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RUNNHP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RUNNOP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RUPPOP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RUPP1P(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RUNNBP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RUNNXP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RUPPBP(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RUPP1B(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WFR(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(CPOOLR(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(ZPOOLR(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(PPOOLR(N,L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTDP1(N,NR,NZ,NY,NX)
     2,NR=1,NRT(NZ,NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTWT1(N,NR,NZ,NY,NX)
     2,NR=1,NRT(NZ,NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTWT1N(N,NR,NZ,NY,NX)
     2,NR=1,NRT(NZ,NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTWT1P(N,NR,NZ,NY,NX)
     2,NR=1,NRT(NZ,NY,NX))
      DO 9965 NR=1,NRT(NZ,NY,NX)
      WRITE(29,94)I,IDATA(3),NZ,(RTN2(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTLG1(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTRT1(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTRT1N(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTRT1P(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTRT2(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTRT2N(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTRT2P(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(RTLG2(N,L,NR,NZ,NY,NX),L=1,NJ(NY,NX))
9965  CONTINUE
9970  CONTINUE
      WRITE(29,94)I,IDATA(3),NZ,(CPOOLN(L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(ZPOOLN(L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(PPOOLN(L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTNDL(L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTNDLN(L,NZ,NY,NX),L=1,NJ(NY,NX))
      WRITE(29,94)I,IDATA(3),NZ,(WTNDLP(L,NZ,NY,NX),L=1,NJ(NY,NX))
      ENDIF
9980  CONTINUE
      ENDIF
9985  CONTINUE
9990  CONTINUE
90    FORMAT(30I4)
91    FORMAT(3I4,4I8,30I4)
92    FORMAT(3I4,10E17.8E3)
93    FORMAT(3I4,15E17.8E3,/,15E17.8E3,/,15E17.8E3,/,15E17.8E3)
94    FORMAT(3I4,21E17.8E3)
95    FORMAT(3I4,10I6)
96    FORMAT(3I4,15E17.8E3,/,15E17.8E3,/,15E17.8E3)
97    FORMAT(3I4,15E17.8E3,/,15E17.8E3,/,15E17.8E3,/,15E17.8E3)
      RETURN
      END

max_stars_repo_name: ThemosTsikas/ECOSYS

max_issues_repo_name: ThemosTsikas/ECOSYS

max_forks_repo_name: ThemosTsikas/ECOSYS


program compute_halo_prop
  use io
  use misc
  use flap, only : command_line_interface
  use compute
  ! use hdf5
  implicit none

  !-------------------------------------
  ! Parameters
  !-------------------------------------
  integer, parameter :: NPART = 2**30
  character(len=200) :: ramses_output_end, ramses_output_start, gal_list_filename, &
       associations_filename, dm_halo_list_filename, brick_file, info_file_end,&
       info_file_start, outfile, halo_to_cpu_file, mergertree_file, param_output_path, &
       param_output
  integer            :: param_from, param_to, param_output_number, param_halo_i
  integer, allocatable :: param_output_number_list(:), param_halos(:)
  integer, parameter :: NPARTICLE_TO_PROBE_HALO = 50, NCPU_PER_HALO = 10
  real(kind=8)       :: param_min_m, param_max_m
  integer            :: param_verbosity
  logical            :: stop_flag

  type(command_line_interface) :: cli
  !-------------------------------------
  ! List data
  !-------------------------------------
  real(kind=4), dimension(:,:), allocatable :: data_gal, data_associations, data_halo
  integer                                   :: ngal, nassoc, ndm_halo, gal_cols, assoc_cols, dm_halo_cols
  !-------------------------------------
  ! Brick data
  !-------------------------------------
  real(kind=8), dimension(:, :), allocatable :: posDM, LDM
  integer, dimension(:), allocatable         :: idDM
  real(kind=8), dimension(:), allocatable    :: rvirDM, mDM, mvirDM, TvirDM, hlevel
  type(INFOS_T)                              :: infos
  type(MEMBERS_T), dimension(:), allocatable :: members
  integer                                    :: nbodies, nb_of_halos, nb_of_subhalos, nDM
  real(kind=4)                               :: aexp_tmp, age_univ
  !-------------------------------------
  ! Merger tree data
  !-------------------------------------
  ! integer,  allocatable, dimension(:)        :: mt_nhalos, mt_nsubhalos
  ! real(kind=4), dimension(:), allocatable    :: mt_aexp, mt_omega_t, mt_age_univ, time
  ! real(kind=4), dimension(:, :), allocatable :: mt_pos, mt_vel, initial_pos
  ! integer                                    :: mt_nsteps, nsteps, nhalos
  ! integer                                    :: initial_halo_i, current_halo, parent_halo, istep, prev, father, step, max_nhalo
  ! integer, dimension(:), allocatable         :: tmp_nhalos, halos_z0
  ! integer, dimension(:), allocatable         :: parent, parent_at_step
  !-------------------------------------
  ! Particle data
  !-------------------------------------
  integer                                    :: ndim, nparts
  real(kind=8), dimension(:, :), allocatable :: vel, tmp_pos
  real(kind=8), dimension(3, NPART)          :: pos, pos_around_halo
  real(kind=4), dimension(:), allocatable    :: birth_date
  real(kind=8), dimension(:), allocatable    :: m
  integer                                    :: nstar, halo_found
  integer,      dimension(:), allocatable    :: tmp_ids
  integer,      dimension(NPART)             :: ids, order, ids_around
  integer, dimension(:), allocatable         :: order2
  !-------------------------------------
  ! Halo to cpu
  !-------------------------------------
  integer                               :: n_cpu_per_halo
  integer, dimension(:, :), allocatable :: halo_to_cpu
  integer, dimension(:), allocatable    :: cpu_list
  !-------------------------------------
  ! Halo properties
  !-------------------------------------
  real(kind=8), dimension(3, 3)               :: I_t_diag
  real(kind=8), dimension(:,:,:), allocatable :: I_t
  real(kind=8), dimension(:), allocatable     :: m_in_halo, m_in_box
  real(kind=8), dimension(:, :), allocatable  :: pos_in_halo, prev_pos_in_halo, vel_in_halo
  integer, dimension(:), allocatable          :: ids_in_box, ids_in_halo, prev_ids_in_halo
  real(kind=8), dimension(:, :), allocatable  :: pos_in_box
  real(kind=8)                                :: mtot
  real(kind=8), dimension(3)                  :: pos_span, pos_mean
  integer                                     :: halo_i, ntot_halo, halo_counter, part_i, parts_in_region
  integer                                     :: part_counter
  integer, dimension(:), allocatable          :: halo_found_mask
  !-------------------------------------
  ! Tmp variables
  !-------------------------------------
  integer                                   :: i, j, k, cpu, i1, i2, counter, x, y, z
  integer                                   :: tmp_int, unit,&
       & out_unit, center_output_unit, tmp_int2, index
  character(len=200)                        :: tmp_char
  real                                      :: tmp_real, factor
  logical                                   :: tmp_bool
  real(kind=8), allocatable, dimension(:)   :: tmp_arr, tmp_arr2
  real(kind=8), allocatable, dimension(:,:) :: tmp_dblarr
  integer, allocatable, dimension(:)        :: tmp_iarr
  logical, dimension(4096)                  :: cpu_read
  real(kind=8), dimension(3)                :: X0, X1, center
  real(kind=8), dimension(:, :), allocatable :: centers
  real(kind=8)                              :: margin
  integer                                   :: nstep, max_nparts, step_i

  !-------------------------------------
  ! random
  !-------------------------------------
  call random_seed()

  !-------------------------------------
  ! Read parameters
  !-------------------------------------
  call cli%init(progname='find_particles_in_halo',&
       description='Find particles around many halos')

  call cli%add(switch='--output-path', help='Path of the output', &
       act='store', def='/data52/Horizon-AGN/OUTPUT_DIR')
  call cli%add(switch='--output-number', help='Number of the output', &
       act='store', def='2', nargs='+')
  call cli%add(switch='--verbose', help='Verbosity', &
       act='store', def='0')
  call cli%add(switch='--halos', help='List of halos to find', act='store', &
       nargs='+', def='1')
  call cli%add(switch='--brick', help='Path for brick file that contains the halos', &
       act='store', def='/data52/Horizon-AGN/TREE_DM_celldx2kpc_SC0.9r/tree_bricks782')
  call cli%add(switch='--output', switch_ab='-o', help='Name of the output file', &
       act='store', def='data.out')

  call cli%get(switch='--output-path', val=param_output_path)
  call cli%get_varying(switch='--output-number', val=param_output_number_list)
  call cli%get(switch='--verbose', val=param_verbosity)
  call cli%get_varying(switch='--halos', val=param_halos)
  call cli%get(switch='--output', val=param_output)

  !-------------------------------------
  ! Read brick file
  !-------------------------------------
  print*, ''
  print*, 'Reading info file'
  call read_info_headers(param_output_path, param_output_number_list(1), infos)

  call cli%get(switch='--brick', val=tmp_char)
  print*, ''
  print*, 'Reading brick file "' // trim(tmp_char) // '"…'
  call read_brick_header(tmp_char, infos, nbodies, aexp_tmp, age_univ,&
       nb_of_halos, nb_of_subhalos)
  nDM = nb_of_halos + nb_of_subhalos
  allocate(idDM(nDM))
  allocate(posDM(infos%ndim, nDM))
  allocate(rvirDM(nDM))
  allocate(mDM(nDM))
  allocate(mvirDM(nDM))
  allocate(TvirDM(nDM))
  allocate(hlevel(nDM))
  allocate(LDM(infos%ndim, nDM))
  allocate(members(nDM))
  allocate(cpu_list(infos%ncpu))
  print*, '    …allocated…'
  call read_brick_data(nDM, infos, .true., &
       & mDM, posDM, rvirDM, mvirDM, TvirDM,&
       & hlevel, LDM, idDM, members)
  print*, '    …red!'
  deallocate(rvirDM, mvirDM, TvirDM, hlevel, LDM)

  counter = 1
  !----------------------------------------
  ! Read all the cpus
  !----------------------------------------
  do cpu = 1, infos%ncpu
     if (param_verbosity >= 3) then
        write(*, '(a, i5, a, i5, a, 1x, i10, a, i10, a)') 'cpu', cpu, '/', infos%ncpu,&
             '(parts', counter, '/', NPART, ')'
     end if

     call read_particle(param_output_path, 2, cpu, nstar, &
          tmp_pos, vel, m, tmp_ids, birth_date, ndim, nparts)
     pos(:, counter:counter + nparts) = tmp_pos
     ids(counter:counter + nparts) = tmp_ids
     ! add the particles into the array
     counter = counter + nparts
  end do

  !----------------------------------------
  ! Sort the indexes
  !----------------------------------------
  if (param_verbosity >= 2) then
     write(*, *) 'Sorting'
  end if
  call quick_sort(ids, order)

  !----------------------------------------
  ! Iterate over each given halo
  !----------------------------------------
  do i = 1, size(param_halos)
     halo_i = param_halos(i)

     if (param_verbosity >= 2) then
        write(*, *) 'Getting halo', halo_i
     end if

     ! allocate data for halo
     allocate(pos_in_halo(infos%ndim, members(halo_i)%parts))

     !----------------------------------------
     ! get the members of the halo
     !----------------------------------------
     do j = 1, members(halo_i)%parts
        tmp_int = indexOf(members(halo_i)%ids(j), ids)

        pos_in_halo(:, j) = pos(:, order(tmp_int))
     end do
     !----------------------------------------
     ! Rectify the position and get center
     !----------------------------------------
     call correct_positions(pos_in_halo)
     call meanval(pos_in_halo, pos_mean, infos%ndim, members(halo_i)%parts)
     pos_span = maxval(pos_in_halo, 2) - minval(pos_in_halo, 2)

     !----------------------------------------
     ! get the particles around halo
     !----------------------------------------
     call filter_region(pos_mean, 2*pos_span, &
          ids, pos, &
          ids_around, pos_around_halo, &
          parts_in_region, order)
     call correct_positions(pos_around_halo(:, :parts_in_region))

     if (param_verbosity >= 3) write(*, *) 'Found', parts_in_region, 'around halo'

     !----------------------------------------
     ! write data
     !----------------------------------------
     write(tmp_char, '(a,i0.5)') param_output_path, halo_i
     if (param_verbosity >= 2) write(*, *) 'Writing output to', trim(tmp_char)
     open(newunit=out_unit, file=trim(tmp_char), form='formatted')

     ! headers
     write(out_unit, '(a9, a5, 6a13)') 'id', 'in', 'x', 'y', 'z'

     ! in halo
     do k = 1, members(halo_i)%parts
        write(out_unit, '(i12, L5, 6ES14.6e2)') &
             members(halo_i)%ids(k), .true., pos_in_halo(:, k)
     end do

     ! out of halo (skipping parts in halo)
     allocate(order2(members(halo_i)%parts))
     call quick_sort(members(halo_i)%ids, order2)

     ! for particles around …
     do k = 1, parts_in_region
        ! … that aren't in the halo …
        if (indexOf(ids_around(k), members(halo_i)%ids) == -1) then
           ! … write their data to the file
           write(out_unit, '(i12, L5, 6ES14.6e2)') &
                ids_around(k), .false., pos_around_halo(:, k)
        end if
     end do
     close(out_unit)

     deallocate(order2, pos_in_halo)

  end do

end program compute_halo_prop


text:
      SUBROUTINE iau_G2ICRS ( DL, DB, DR, DD )
*+
*  - - - - - - - - - - -
*   i a u _ G 2 I C R S
*  - - - - - - - - - - -
*
*  Transformation from Galactic Coordinates to ICRS.
*
*  This routine is part of the International Astronomical Union's
*  SOFA (Standards of Fundamental Astronomy) software collection.
*
*  Status:  support routine.
*
*  Given:
*     DL       d      galactic longitude (radians)
*     DB       d      galactic latitude (radians)
*
*  Returned:
*     DR       d      ICRS right ascension (radians)
*     DD       d      ICRS declination (radians)
*
*  Notes:
*
*  1) The IAU 1958 system of Galactic coordinates was defined with
*     respect to the now obsolete reference system FK4 B1950.0.  When
*     interpreting the system in a modern context, several factors have
*     to be taken into account:
*
*     . The inclusion in FK4 positions of the E-terms of aberration.
*
*     . The distortion of the FK4 proper motion system by differential
*       Galactic rotation.
*
*     . The use of the B1950.0 equinox rather than the now-standard
*       J2000.0.
*
*     . The frame bias between ICRS and the J2000.0 mean place system.
*
*     The Hipparcos Catalogue (Perryman & ESA 1997) provides a rotation
*     matrix that transforms directly between ICRS and Galactic
*     coordinates with the above factors taken into account.  The
*     matrix is derived from three angles, namely the ICRS coordinates
*     of the Galactic pole and the longitude of the ascending node of
*     the galactic equator on the ICRS equator.  They are given in
*     degrees to five decimal places and for canonical purposes are
*     regarded as exact.  In the Hipparcos Catalogue the matrix elements
*     are given to 10 decimal places (about 20 microarcsec).  In the
*     present SOFA routine the matrix elements have been recomputed from
*     the canonical three angles and are given to 30 decimal places.
*
*  2) The inverse transformation is performed by the routine iau_ICRS2G.
*
*  Called:
*     iau_ANP      normalize angle into range 0 to 2pi
*     iau_ANPM     normalize angle into range +/- pi
*     iau_S2C      spherical coordinates to unit vector
*     iau_TRXP     product of transpose of r-matrix and p-vector
*     iau_C2S      p-vector to spherical
*
*  Reference:
*     Perryman M.A.C. & ESA, 1997, ESA SP-1200, The Hipparcos and Tycho
*     catalogues.  Astrometric and photometric star catalogues
*     derived from the ESA Hipparcos Space Astrometry Mission.  ESA
*     Publications Division, Noordwijk, Netherlands.
*
*  This revision:   2015 January 9
*
*  SOFA release 2020-07-21
*
*  Copyright (C) 2020 IAU SOFA Board.  See notes at end.
*
*-----------------------------------------------------------------------

      IMPLICIT NONE
      DOUBLE PRECISION DL, DB, DR, DD

      DOUBLE PRECISION iau_ANP, iau_ANPM

      DOUBLE PRECISION V1(3), V2(3)

*
*  L2,B2 system of galactic coordinates in the form presented in the
*  Hipparcos Catalogue.  In degrees:
*
*  P = 192.85948    right ascension of the Galactic north pole in ICRS
*  Q =  27.12825    declination of the Galactic north pole in ICRS
*  R =  32.93192    longitude of the ascending node of the Galactic
*                   plane on the ICRS equator
*
*  ICRS to galactic rotation matrix, obtained by computing
*  R_3(-R) R_1(pi/2-Q) R_3(pi/2+P) to the full precision shown:
*
      DOUBLE PRECISION R(3,3)
      DATA R(1,1), R(1,2), R(1,3),
     :     R(2,1), R(2,2), R(2,3),
     :     R(3,1), R(3,2), R(3,3) /
     :    -0.054875560416215368492398900454D0,
     :    -0.873437090234885048760383168409D0,
     :    -0.483835015548713226831774175116D0,
     :    +0.494109427875583673525222371358D0,
     :    -0.444829629960011178146614061616D0,
     :    +0.746982244497218890527388004556D0,
     :    -0.867666149019004701181616534570D0,
     :    -0.198076373431201528180486091412D0,
     :    +0.455983776175066922272100478348D0 /

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*  Spherical to Cartesian.
      CALL iau_S2C ( DL, DB, V1 )

*  Galactic to ICRS.
      CALL iau_TRXP ( R, V1, V2 )

*  Cartesian to spherical.
      CALL iau_C2S ( V2, DR, DD )

*  Express in conventional ranges.
      DR = iau_ANP ( DR )
      DD = iau_ANPM ( DD )

*  Finished.

*+----------------------------------------------------------------------
*
*  Copyright (C) 2020
*  Standards Of Fundamental Astronomy Board
*  of the International Astronomical Union.
*
*  =====================
*  SOFA Software License
*  =====================
*
*  NOTICE TO USER:
*
*  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
*  CONDITIONS WHICH APPLY TO ITS USE.
*
*  1. The Software is owned by the IAU SOFA Board ("SOFA").
*
*  2. Permission is granted to anyone to use the SOFA software for any
*     purpose, including commercial applications, free of charge and
*     without payment of royalties, subject to the conditions and
*     restrictions listed below.
*
*  3. You (the user) may copy and distribute SOFA source code to others,
*     and use and adapt its code and algorithms in your own software,
*     on a world-wide, royalty-free basis.  That portion of your
*     distribution that does not consist of intact and unchanged copies
*     of SOFA source code files is a "derived work" that must comply
*     with the following requirements:
*
*     a) Your work shall be marked or carry a statement that it
*        (i) uses routines and computations derived by you from
*        software provided by SOFA under license to you; and
*        (ii) does not itself constitute software provided by and/or
*        endorsed by SOFA.
*
*     b) The source code of your derived work must contain descriptions
*        of how the derived work is based upon, contains and/or differs
*        from the original SOFA software.
*
*     c) The names of all routines in your derived work shall not
*        include the prefix "iau" or "sofa" or trivial modifications
*        thereof such as changes of case.
*
*     d) The origin of the SOFA components of your derived work must
*        not be misrepresented;  you must not claim that you wrote the
*        original software, nor file a patent application for SOFA
*        software or algorithms embedded in the SOFA software.
*
*     e) These requirements must be reproduced intact in any source
*        distribution and shall apply to anyone to whom you have
*        granted a further right to modify the source code of your
*        derived work.
*
*     Note that, as originally distributed, the SOFA software is
*     intended to be a definitive implementation of the IAU standards,
*     and consequently third-party modifications are discouraged.  All
*     variations, no matter how minor, must be explicitly marked as
*     such, as explained above.
*
*  4. You shall not cause the SOFA software to be brought into
*     disrepute, either by misuse, or use for inappropriate tasks, or
*     by inappropriate modification.
*
*  5. The SOFA software is provided "as is" and SOFA makes no warranty
*     as to its use or performance.   SOFA does not and cannot warrant
*     the performance or results which the user may obtain by using the
*     SOFA software.  SOFA makes no warranties, express or implied, as
*     to non-infringement of third party rights, merchantability, or
*     fitness for any particular purpose.  In no event will SOFA be
*     liable to the user for any consequential, incidental, or special
*     damages, including any lost profits or lost savings, even if a
*     SOFA representative has been advised of such damages, or for any
*     claim by any third party.
*
*  6. The provision of any version of the SOFA software under the terms
*     and conditions specified herein does not imply that future
*     versions will also be made available under the same terms and
*     conditions.
*
*  In any published work or commercial product which uses the SOFA
*  software directly, acknowledgement (see www.iausofa.org) is
*  appreciated.
*
*  Correspondence concerning SOFA software should be addressed as
*  follows:
*
*      By email:  sofa@ukho.gov.uk
*      By post:   IAU SOFA Center
*                 HM Nautical Almanac Office
*                 UK Hydrographic Office
*                 Admiralty Way, Taunton
*                 Somerset, TA1 2DN
*                 United Kingdom
*
*-----------------------------------------------------------------------

      END

ext: for

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
      LOGICAL FUNCTION AB13ID( JOBSYS, JOBEIG, EQUIL, CKSING, RESTOR,
     $                         UPDATE, N, M, P, A, LDA, E, LDE, B, LDB,
     $                         C, LDC, NR, RANKE, TOL, IWORK, DWORK,
     $                         LDWORK, IWARN, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To check whether the transfer function
C
C                                     -1
C       G(lambda) := C*( lambda*E - A ) *B
C
C     of a given linear time-invariant descriptor system with
C     generalized state space realization (lambda*E-A,B,C) is proper.
C     Optionally, if JOBEIG = 'A', the system (lambda*E-A,B,C) is
C     reduced to an equivalent one (lambda*Er-Ar,Br,Cr) with only
C     controllable and observable eigenvalues in order to use it for a
C     subsequent L_inf-norm computation; if JOBEIG = 'I', the system is
C     reduced to an equivalent one (lambda*Er-Ar,Br,Cr) without
C     uncontrollable and unobservable infinite eigenvalues. In this
C     case, intended mainly for checking the properness, the returned
C     system is not fully reduced, unless UPDATE = 'U'.
C
C     FUNCTION VALUE
C
C     AB13ID  LOGICAL
C             Indicates whether the transfer function is proper.
C             If AB13ID = .TRUE., the transfer function is proper;
C             otherwise, it is improper.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOBSYS  CHARACTER*1
C             Indicates whether the system (lambda*E-A,B,C) is already
C             in the reduced form which is obtained as stated in
C             JOBEIG, as follows.
C             = 'R': The system is not in a reduced form, the reduction
C                    step is performed;
C             = 'N': The system is in a reduced form; the reduction step
C                    is omitted.
C
C     JOBEIG  CHARACTER*1
C             Indicates which kind of eigenvalues of the matrix pencil
C             lambda*E-A should be removed if JOBSYS = 'R', as follows:
C             = 'A': All uncontrollable and unobservable eigenvalues
C                    are removed; the reduced system is returned in
C                    the arrays A, E, B, C;
C             = 'I': Only all uncontrollable and unobservable infinite
C                    eigenvalues are removed; the returned system is not
C                    fully reduced if UPDATE = 'N'.
C
C     EQUIL   CHARACTER*1
C             Specifies whether the user wishes to preliminarily scale
C             the system (lambda*E-A,B,C) as follows:
C             = 'S': Perform scaling;
C             = 'N': Do not perform scaling.
C
C     CKSING  CHARACTER*1
C             Specifies whether the user wishes to check if the pencil
C             (lambda*E-A) is singular as follows:
C             = 'C':  Check singularity;
C             = 'N':  Do not check singularity.
C             If the pencil is singular, the reduced system computed for
C             CKSING = 'N' may have completely different eigenvalues
C             than the given system.
C             The test is performed only if JOBSYS = 'R'.
C
C     RESTOR  CHARACTER*1
C             Specifies whether the user wishes to save the system
C             matrices before each reduction phase (if JOBSYS = 'R') and
C             restore them if no order reduction took place as follows:
C             = 'R':  Save and restore;
C             = 'N':  Do not save the matrices.
C             This option is ineffective if JOBSYS = 'N'.
C
C     UPDATE  CHARACTER*1
C             Specifies whether the user wishes to update the matrices
C             A, B, and C if JOBEIG = 'I' as follows:
C             = 'U':  Update the matrices A, B and C;
C             = 'N':  Do not update the matrices A, B and C when
C                     performing URV decomposition of the matrix E
C                     (see METHOD).
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The dimension of the descriptor state vector; also the
C             order of square matrices A and E, the number of rows of
C             matrix B, and the number of columns of matrix C.  N >= 0.
C
C     M       (input) INTEGER
C             The dimension of the descriptor system input vector; also
C             the number of columns of matrix B.  M >= 0.
C
C     P       (input) INTEGER
C             The dimension of the descriptor system output vector; also
C             the number of rows of matrix C.  P >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array must
C             contain the state matrix A.
C             On exit, if JOBSYS = 'R' and JOBEIG = 'A', the leading
C             NR-by-NR part of this array contains the reduced order
C             state matrix Ar of a fully controllable and observable
C             realization for the original system. If JOBSYS = 'R' and
C             JOBEIG = 'I', the leading NR-by-NR part of this array
C             contains the reduced order state matrix Ar of a
C             transformed system without uncontrollable and unobservable
C             infinite poles. In this case, the matrix Ar does not
C             correspond to the returned matrix Er (obtained after a
C             URV decomposition), unless UPDATE = 'U' or RANKE < NR.
C             On exit, if JOBSYS = 'N' and (JOBEIG = 'A' or UPDATE = 'U'
C             or RANKE < N), the leading N-by-N part of this array
C             contains the transformed matrix A corresponding to the
C             URV decomposition of E (see (2) in METHOD), and if
C             JOBEIG = 'I' and UPDATE = 'N', the submatrix A22 in (2) is
C             further transformed to estimate its rank.
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1,N).
C
C     E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
C             On entry, the leading N-by-N part of this array must
C             contain the descriptor matrix E.
C             On exit, if JOBSYS = 'R' and JOBEIG = 'A', the leading
C             NR-by-NR part of this array contains the reduced order
C             descriptor matrix Er of a completely controllable and
C             observable realization for the original system. The
C             reduced matrix Er is in upper triangular form.
C             If JOBSYS = 'R' and JOBEIG = 'I', the leading NR-by-NR
C             part of this array contains the reduced order descriptor
C             matrix Er of a transformed system without uncontrollable
C             and unobservable infinite poles. The reduced matrix Er is
C             upper triangular. In both cases, or if JOBSYS = 'N', the
C             matrix Er results from a URV decomposition of the matrix E
C             (see METHOD).
C
C     LDE     INTEGER
C             The leading dimension of the array E.  LDE >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension
C             (LDB,MAX(M,P))
C             On entry, the leading N-by-M part of this array must
C             contain the input matrix B; the remainder of the leading
C             N-by-MAX(M,P) part is used as internal workspace.
C             On exit, if JOBSYS = 'R' and JOBEIG = 'A', the leading
C             NR-by-M part of this array contains the reduced input
C             matrix Br of a completely controllable and observable
C             realization for the original system. If JOBSYS = 'R' and
C             JOBEIG = 'I', the leading NR-by-M part of this array
C             contains the transformed input matrix Br obtained after
C             removing the uncontrollable and unobservable infinite
C             poles; the transformations for the URV decomposition of
C             the matrix E are not applied if UPDATE = 'N'.
C             On exit, if JOBSYS = 'N' and (JOBEIG = 'A' or
C             UPDATE = 'U'), the leading N-by-M part of this array
C             contains the transformed matrix B corresponding to the
C             URV decomposition of E, but if JOBEIG = 'I', EQUIL = 'N'
C             and UPDATE = 'N', the array B is unchanged on exit.
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= MAX(1,N).
C
C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
C             On entry, the leading P-by-N part of this array must
C             contain the output matrix C; the remainder of the leading
C             MAX(M,P)-by-N part is used as internal workspace.
C             On exit, if JOBSYS = 'R' and JOBEIG = 'A', the leading
C             P-by-NR part of this array contains the transformed output
C             matrix Cr of a completely controllable and observable
C             realization for the original system. If JOBSYS = 'R' and
C             JOBEIG = 'I', the leading P-by-NR part of this array
C             contains the transformed output matrix Cr obtained after
C             removing the uncontrollable and unobservable infinite
C             poles; the transformations for the URV decomposition of
C             the matrix E are not applied if UPDATE = 'N'.
C             On exit, if JOBSYS = 'N' and (JOBEIG = 'A' or
C             UPDATE = 'U'), the leading P-by-N part of this array
C             contains the transformed matrix C corresponding to the
C             URV decomposition of E, but if JOBEIG = 'I', EQUIL = 'N'
C             and UPDATE = 'N', the array C is unchanged on exit.
C
C     LDC     INTEGER
C             The leading dimension of the array C.
C             LDC >= MAX(1,M,P) if N > 0;
C             LDC >= 1          if N = 0.
C
C     NR      (output) INTEGER
C             The order of the reduced generalized state space
C             representation (lambda*Er-Ar,Br,Cr) as stated in JOBEIG.
C             If JOBEIG = 'A', NR denotes the order of a reduced system
C             without any uncontrollable or unobservable eigenvalues; if
C             JOBEIG = 'I', NR denotes the order of the reduced system
C             without any uncontrollable or unobservable infinite
C             eigenvalues. If JOBSYS = 'N', then NR = N.
C
C     RANKE   (output) INTEGER
C             The effective (estimated) rank of the reduced matrix Er.
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION array, dimension 3
C             TOL(1) is the tolerance to be used in rank determinations
C             when transforming (lambda*E-A,B,C). If the user sets
C             TOL(1) > 0, then the given value of TOL(1) is used as a
C             lower bound for reciprocal condition numbers in rank
C             determinations; a (sub)matrix whose estimated condition
C             number is less than 1/TOL(1) is considered to be of full
C             rank.  If the user sets TOL(1) <= 0, then an implicitly
C             computed, default tolerance, defined by TOLDEF1 = N*N*EPS,
C             is used instead, where EPS is the machine precision (see
C             LAPACK Library routine DLAMCH).  TOL(1) < 1.
C             TOL(2) is the tolerance to be used for checking pencil
C             singularity when CKSING = 'C', or singularity of the
C             matrices A and E when CKSING = 'N'. If the user sets
C             TOL(2) > 0, then the given value of TOL(2) is used.
C             If the user sets TOL(2) <= 0, then an implicitly
C             computed, default tolerance, defined by  TOLDEF2 = 10*EPS,
C             is used instead.  TOL(2) < 1.
C             TOL(3) is the threshold value for magnitude of the matrix
C             elements, if EQUIL = 'S': elements with magnitude less
C             than or equal to TOL(3) are ignored for scaling. If the
C             user sets TOL(3) >= 0, then the given value of TOL(3) is
C             used. If the user sets TOL(3) < 0, then an implicitly
C             computed, default threshold, defined by  THRESH = c*EPS,
C             where c = MAX(norm_1(A,E,B,C)) is used instead.
C             TOL(3) = 0 is not always a good choice.  TOL(3) < 1.
C             TOL(3) is not used if EQUIL = 'N'.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension LIWORK
C             If JOBSYS = 'R',  LIWORK >= 2*N+MAX(M,P)+7;
C             If JOBSYS = 'N',  LIWORK >= N.
C             If JOBSYS = 'R', the first 7 elements of IWORK contain
C             information on performed reduction and on structure of
C             resulting system matrices after removing the specified
C             eigenvalues (see the description of the parameter INFRED
C             of the SLICOT Library routine TG01JY).
C
C     DWORK   DOUBLE PRECISION array, dimension LDWORK
C             On exit, if  INFO = 0,  DWORK(1) returns the optimal value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             If JOBSYS = 'R', and EQUIL = 'S',
C                LDWORK >= MAX(w+4*N+4,8*N,x,y),
C             where w = N*N,                  if JOBEIG = 'A',
C                   w = 0,                    if JOBEIG = 'I',
C                   x = MAX(2*(z+MAX(M,P)+N-1),N*N+4*N), if RESTOR = 'R'
C                   x = MAX(  2*(MAX(M,P)+N-1),N*N+4*N), if RESTOR = 'N'
C                   y = 2*N*N+10*N+MAX(N,23), if CKSING = 'C',
C                   y = 0,                    if CKSING = 'N',
C                   z = 2*N*N+N*M+N*P;
C             if JOBSYS = 'R', and EQUIL = 'N',
C                LDWORK >= MAX(w+4*N+4,x,y);
C             if JOBSYS = 'N', and JOBEIG = 'A'  or UPDATE = 'U',
C                              and EQUIL  = 'S',
C                LDWORK >= MAX(N*N+4*N+4,8*N,N+M,N+P);
C             if JOBSYS = 'N', and JOBEIG = 'A'  or UPDATE = 'U',
C                              and EQUIL  = 'N',
C                LDWORK >= MAX(N*N+4*N+4,N+M,N+P);
C             if JOBSYS = 'N', and JOBEIG = 'I' and UPDATE = 'N',
C                              and EQUIL  = 'S',
C                LDWORK >= MAX(4*N+4,8*N);
C             if JOBSYS = 'N', and JOBEIG = 'I' and UPDATE = 'N',
C                              and EQUIL  = 'N',
C                LDWORK >= 4*N+4.
C             If JOBSYS = 'R' and ( RESTOR = 'R' or
C             LDWORK >= MAX(1,2*N*N+N*M+N*P+2*(MAX(M,P)+N-1) ),
C             then more accurate results are to be expected by
C             considering only those reduction phases in the SLICOT
C             Library routine TG01JY, where effective order reduction
C             occurs. This is achieved by saving the system matrices
C             before each phase (after orthogonally triangularizing the
C             matrix A or the matrix E, if RESTOR = 'N') and restoring
C             them if no order reduction took place. However, higher
C             global accuracy is not guaranteed.
C             For good performance, LDWORK should be generally larger.
C
C             If LDWORK = -1, then a workspace query is assumed;
C             the routine only calculates the optimal size of the
C             DWORK array, returns this value as the first entry of
C             the DWORK array, and no error message related to LDWORK
C             is issued by XERBLA. The optimal workspace includes the
C             extra space for improving the accuracy.
C
C     Error/Warning Indicator
C
C     IWARN   INTEGER
C             = 0: When determining the rank of a matrix, the distance
C                  between the tolerance TOL(1) and the estimated
C                  singular values is sufficiently large. The rank can
C                  be safely determined;
C             = 1: When determining the rank of a matrix, there exist
C                  estimated singular values which are very close to the
C                  tolerance TOL(1). The computed rank is possibly
C                  incorrect.
C
C     INFO    INTEGER
C             = 0: succesful exit;
C             < 0: if INFO = -i, the i-th argument had an illegal value;
C             = 1: the given pencil A - lambda*E is numerically
C                  singular and the reduced system is not computed.
C                  However, the system is considered improper, and
C                  AB13ID is set to .FALSE.
C                  This error can be returned only if CKSING = 'C'.
C
C     METHOD
C
C     If JOBSYS = 'R', the routine first removes uncontrollable and
C     unobservable infinite eigenvalues of the pencil lambda*E-A. If, in
C     addition, JOBEIG = 'A', uncontrollable and unobservable zero
C     eigenvalues are also removed. Then, or if JOBSYS = 'N', a
C     URV decomposition of the matrix E is performed, i.e., orthogonal
C     matrices U and V are computed, such that
C
C               ( T  0 )
C       U*E*V = (      ) with a full-rank matrix T.                  (1)
C               ( 0  0 )
C
C     Then the matrix A (or a copy of A if JOBEIG = 'A' or UPDATE = 'U')
C     is updated and partioned as in (1), i.e.,
C
C               ( A11  A12 )
C       U*A*V = (          ) ,                                       (2)
C               ( A21  A22 )
C
C     and the rank of A22 is computed. If A22 is invertible, the
C     transfer function is proper, otherwise it is improper. If
C     required (i.e., JOBEIG = 'A' or UPDATE = 'U'), the matrices B and
C     C are updated as well in order to obtain an equivalent reduced
C     system with the same transfer function. See also Chapter 3 in [1],
C     [2] for more details.
C
C     REFERENCES
C
C     [1] Voigt, M.
C         L_inf-Norm Computation for Descriptor Systems.
C         Diploma Thesis, Chemnitz University of Technology, Department
C         of Mathematics, Germany, July 2010.
C
C     [2] Benner, P., Sima, V., Voigt, M.
C         L_infinity-norm computation for continuous-time descriptor
C         systems using structured matrix pencils.
C         IEEE Trans. Automat. Contr., vol. 57, pp. 233-238, 2012.
C
C     NUMERICAL ASPECTS
C
C     The algorithm requires O(N**3) floating point operations. During
C     the algorithm it is necessary to determine the rank of certain
C     matrices. Therefore it is crucial to use an appropriate tolerance
C     TOL(1) to make correct rank decisions.
C
C     CONTRIBUTORS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Romania,
C     Jan. 2012.
C     Based on the subroutine DGEISP by Matthias Voigt, Chemnitz
C     University of Technology, Department of Mathematics, Feb. 2010.
C
C     REVISIONS
C
C     V. Sima, Feb. 2012, March 2012, April 2012, June 2012.
C
C     KEYWORDS
C
C     Descriptor system, proper transfer function, L_inf-norm
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TEN
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1 )
C
C     .. Scalar Arguments ..
      CHARACTER          CKSING, EQUIL, JOBEIG, JOBSYS, RESTOR, UPDATE
      INTEGER            INFO, IWARN, LDA, LDB, LDC, LDE, LDWORK, M, N,
     $                   NR, P, RANKE
C
C     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   DWORK( * ), E( LDE, * ), TOL( * )
C
C     .. Local Scalars ..
      CHARACTER          SYSTYP
      LOGICAL            LEQUIL, LQUERY, LREDC, LREMA, LRUPD, LSING,
     $                   LUPD, MAXACC
      INTEGER            I, ISV, ITAU, IWRK, IWS, J, K, MAXMP, MAXWRK,
     $                   MINWRK, N1, NA, RANKA
      DOUBLE PRECISION   PREC, SVLMAX, THRESH, TOLDEF
C
C     .. Local Arrays ..
      DOUBLE PRECISION   DUM( 2 ), TOLV( 3 )
C
C     .. External Subroutines ..
      EXTERNAL           DLACPY, DLASET, DORMQR, DORMRZ, DSWAP, DTZRZF,
     $                   MB03OD, TG01AD, TG01JY, XERBLA
C
C     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           DLAMCH, DLANGE, LSAME
C
C     .. Intrinsic Functions ..
      INTRINSIC          ABS, INT, MAX, MIN
C
C     .. Executable Statements ..
C
      INFO  = 0
      IWARN = 0
      MAXMP = MAX( M, P )
      N1    = MAX( 1, N )
C
C     Decode the input arguments.
C
      LREDC  = LSAME( JOBSYS, 'R' )
      LREMA  = LSAME( JOBEIG, 'A' )
      LEQUIL = LSAME( EQUIL,  'S' )
      LSING  = LSAME( CKSING, 'C' )
      MAXACC = LSAME( RESTOR, 'R' )
      LUPD   = LSAME( UPDATE, 'U' )
      LRUPD  = LREMA .OR. LUPD
C
C     Test the input arguments.
C
      IF( .NOT.LREDC .AND. .NOT.LSAME( JOBSYS, 'N' ) ) THEN
         INFO = -1
      ELSE IF( .NOT.LREMA  .AND. .NOT.LSAME( JOBEIG, 'I' ) ) THEN
         INFO = -2
      ELSE IF( .NOT.LEQUIL .AND. .NOT.LSAME( EQUIL,  'N' ) ) THEN
         INFO = -3
      ELSE IF( .NOT.LSING  .AND. .NOT.LSAME( CKSING, 'N' ) ) THEN
         INFO = -4
      ELSE IF( .NOT.MAXACC .AND. .NOT.LSAME( RESTOR, 'N' ) ) THEN
         INFO = -5
      ELSE IF( .NOT.LUPD   .AND. .NOT.LSAME( UPDATE, 'N' ) ) THEN
         INFO = -6
      ELSE IF( N.LT.0 ) THEN
         INFO = -7
      ELSE IF( M.LT.0 ) THEN
         INFO = -8
      ELSE IF( P.LT.0 ) THEN
         INFO = -9
      ELSE IF( LDA.LT.N1 ) THEN
         INFO = -11
      ELSE IF( LDE.LT.N1 ) THEN
         INFO = -13
      ELSE IF( LDB.LT.N1 ) THEN
         INFO = -15
      ELSE IF( LDC.LT.1 .OR. ( N.GT.0 .AND. LDC.LT.MAXMP ) ) THEN
         INFO = -17
      ELSE IF( TOL( 1 ).GE.ONE ) THEN
         INFO = -20
      ELSE IF( TOL( 2 ).GE.ONE ) THEN
         INFO = -20
      ELSE IF( LEQUIL ) THEN
         THRESH = TOL( 3 )
         IF( THRESH.GE.ONE )
     $      INFO = -20
      END IF
      IF( INFO.EQ.0 ) THEN
C
C        Compute minimal workspace.
C
         IF( LREDC ) THEN
            K = N*( 2*N + M + P )
            IF( MAXACC ) THEN
               MINWRK = MAX( 1, 2*( K + MAXMP + N - 1 ) )
            ELSE
               MINWRK = MAX( 1, 2*( MAXMP + N - 1 ) )
            END IF
            MINWRK = MAX( MINWRK, N*N + 4*N )
            IF( LSING )
     $         MINWRK = MAX( MINWRK, 2*N*N + 10*N + MAX( N, 23 ) )
C
            IF( LREMA ) THEN
               SYSTYP = 'R'
            ELSE
               SYSTYP = 'P'
            END IF
         ELSE
            MINWRK = 0
         END IF
         IF( LRUPD ) THEN
            MINWRK = MAX( MINWRK, N*N + 4*N + 4, N + MAXMP )
         ELSE
            MINWRK = MAX( MINWRK, 4*N + 4 )
         END IF
         IF( LEQUIL )
     $      MINWRK = MAX( MINWRK, 8*N )
C
         MAXWRK = MINWRK
         LQUERY = LDWORK.EQ.-1
         IF( LQUERY ) THEN
C
C           Compute optimal workspace.
C
            IF( LREDC ) THEN
               CALL TG01JY( 'Irreducible', SYSTYP, 'No Scaling', CKSING,
     $                      RESTOR, N, M, P, A, LDA, E, LDE, B, LDB, C,
     $                      LDC, NR, IWORK, TOL, IWORK, DWORK, -1,
     $                      INFO )
               MAXWRK = MAX( MAXWRK, INT( DWORK( 1 ) ) )
            END IF
            CALL MB03OD( 'QR Decomposition', N, N, E, LDE, IWORK, TOL,
     $                   ZERO, DWORK, RANKE, DWORK, DUM, -1, INFO )
            MAXWRK = MAX( MAXWRK, INT( DUM( 1 ) ) + N + 3 )
            CALL DORMQR( 'Left', 'Transpose', N, N, N, E, LDE, DWORK, A,
     $                   LDA, DUM, -1, INFO )
            MAXWRK = MAX( MAXWRK, INT( DUM( 1 ) ) + N )
            IF( LRUPD ) THEN
               CALL DORMQR( 'Left', 'Transpose', N, M, N, E, LDE, DWORK,
     $                      B, LDB, DUM( 2 ), -1, INFO )
               MAXWRK = MAX( MAXWRK, INT( DUM( 2 ) ) + N )
            END IF
            CALL DTZRZF( N, N, E, LDE, DWORK, DUM, -1, INFO )
            CALL DORMRZ( 'Right', 'Transpose', N, N, N, N, E, LDE,
     $                   DWORK, A, LDA, DUM( 2 ), -1, INFO )
            MAXWRK = MAX( MAXWRK, MAX( INT( DUM( 1 ) ),
     $                                 INT( DUM( 2 ) ) ) + N )
            IF( LRUPD ) THEN
               CALL DORMRZ( 'Right', 'Transpose', P, N, N, N, E, LDE,
     $                      DWORK, C, LDC, DUM, -1, INFO )
               MAXWRK = MAX( MAXWRK, INT( DUM( 1 ) ) + N )
            END IF
         ELSE IF( LDWORK.LT.MINWRK ) THEN
            INFO = -23
         END IF
      END IF
C
      IF( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'AB13ID', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         DWORK( 1 ) = MAXWRK
         RETURN
      END IF
C
C     Quick return if possible.
C
      NR = N
      IF( N.EQ.0 ) THEN
         RANKE  = 0
         AB13ID = .TRUE.
         DWORK( 1 ) = ONE
         RETURN
      END IF
C
C     Set the tolerances.
C
      TOLDEF = TOL( 1 )
      IF( TOLDEF.LE.ZERO .OR. LEQUIL ) THEN
         PREC = DLAMCH( 'Precision' )
         IF( LEQUIL ) THEN
            IF( THRESH.LT.ZERO )
     $         THRESH = MAX( DLANGE( '1-norm', N, N, A, LDA, DWORK ),
     $                       DLANGE( '1-norm', N, N, E, LDE, DWORK ),
     $                       DLANGE( '1-norm', N, M, B, LDB, DWORK ),
     $                       DLANGE( '1-norm', P, N, C, LDC, DWORK ) )*
     $                  PREC
         ELSE
            TOLDEF = N*N*PREC
         END IF
      END IF
      TOLV( 1 ) = TOLDEF
      TOLV( 2 ) = TOL( 2 )
      TOLV( 3 ) = TOL( 3 )
C
C     Computations.
C
C     (Note: Comments in the code beginning "Workspace:" describe the
C     minimal amount of real workspace needed at that point in the
C     code, as well as the preferred amount for good performance.
C     NB refers to the optimal block size for the immediately
C     following subroutine, as returned by ILAENV.)
C
C     Equilibrate the system, if required.
C
C     Workspace: need   8*N.
C
      IF( LEQUIL ) THEN
         CALL TG01AD( 'All', N, N, M, P, THRESH, A, LDA, E, LDE, B, LDB,
     $                C, LDC, DWORK, DWORK( N+1 ), DWORK( 2*N+1 ),
     $                INFO )
         MAXWRK = MAX( MAXWRK, 8*N )
      END IF
C
      IF( LREDC ) THEN
C
C        Step 1: If JOBEIG = 'A' remove all uncontrollable and
C                unobservable eigenvalues of the system; otherwise
C                remove only the uncontrollable and unobservable
C                infinite eigenvalues explicitely, uncontrollable or
C                unobservable zero eigenvalues are not removed.
C
C        Workspace: need   MAX(2*(MAX(M,P)+N-1),N*N+4*N,y),
C                   prefer MAX(2*N*N+N*M+N*P+2*(MAX(M,P)+N-1),y), or
C                          even larger (see, optimal space above).
C
         IWS = 8
         CALL TG01JY( 'Irreducible', SYSTYP, 'No Scaling', CKSING,
     $                RESTOR, N, M, P, A, LDA, E, LDE, B, LDB, C, LDC,
     $                NR, IWORK, TOLV, IWORK( IWS ), DWORK, LDWORK,
     $                INFO )
         MAXWRK = MAX( MAXWRK, INT( DWORK( 1 ) ) )
         IF( INFO.EQ.1 ) THEN
            AB13ID = .FALSE.
            RETURN
         END IF
      ELSE
         IWS = 1
      END IF
C
C     Step 2: Compute a URV decomposition of E(1:NR,1:NR).
C
      DO 10 I = IWS, IWS + NR - 1
         IWORK( I ) = 0
   10 CONTINUE
C
C     Perform a rank-revealing QR decomposition of E(1:NR,1:NR) with
C     respect to the tolerance TOLDEF.
C
C     Workspace: need   4*NR+4,
C                prefer 3*NR+( NR+1 )*NB+3, with NB for DGEQP3.
C
      SVLMAX = DLANGE( 'Frobenius', N, N, E, LDE, DWORK )
      CALL MB03OD( 'QR Decomposition', NR, NR, E, LDE, IWORK( IWS ),
     $             TOLDEF, SVLMAX, DWORK, RANKE, DWORK( NR+1 ),
     $             DWORK( NR+4 ), LDWORK-( NR+3 ), INFO )
      MAXWRK = MAX( MAXWRK, INT( DWORK( NR+4 ) ) + NR + 3 )
C
C     Check whether there are singular values close to the tolerance
C     TOLDEF.
C
      IF( MIN( ABS( DWORK( NR+2 ) - TOLDEF ),
     $         ABS( DWORK( NR+3 ) - TOLDEF ) ) .LT.TOLDEF/TEN )
C
C        One singular value is close to TOLDEF, warning returned.
C
     $   IWARN = 1
C
      IF( RANKE.LT.NR .OR. LRUPD ) THEN
C
C        Perform the same transformations on A(1:NR,1:NR).
C
C        Workspace: need   2*NR,
C                   prefer NR+NR*( NB+1 ).
C
         CALL DORMQR( 'Left', 'Transpose', NR, NR, NR, E, LDE, DWORK, A,
     $                LDA, DWORK( NR+1 ), LDWORK-NR, INFO )
         MAXWRK = MAX( MAXWRK, INT( DWORK( NR+1 ) ) + NR )
C
         DO 20 I = IWS, IWS + NR - 1
            IWORK( I ) = -IWORK( I )
   20    CONTINUE
C
         DO 40 I = IWS, IWS + NR - 1
            IF( IWORK( I ).LT.0 ) THEN
               J = I
               IWORK( J ) = -IWORK( J )
   30          CONTINUE
               K = IWORK( J ) + IWS - 1
               IF( IWORK( K ).LT.0 ) THEN
                  CALL DSWAP( NR, A( 1, J-IWS+1 ), 1, A( 1, K-IWS+1 ),
     $                        1 )
                  IWORK( K ) = -IWORK( K )
                  J = K
                  GO TO 30
               END IF
            END IF
   40    CONTINUE
C
         IF( LRUPD ) THEN
C
C           Update B and C.
C
C           Workspace: need   NR+M,
C                      prefer NR+M*NB.
C
            CALL DORMQR( 'Left', 'Transpose', NR, M, NR, E, LDE, DWORK,
     $                   B, LDB, DWORK( NR+1 ), LDWORK-NR, INFO )
            MAXWRK = MAX( MAXWRK, INT( DWORK( NR+1 ) ) + NR )
C
            DO 50 I = IWS, IWS + NR - 1
               IWORK( I ) = -IWORK( I )
   50       CONTINUE
C
            DO 70 I = IWS, IWS + NR - 1
               IF( IWORK( I ).LT.0 ) THEN
                  J = I
                  IWORK( J ) = -IWORK( J )
   60             CONTINUE
                  K = IWORK( J ) + IWS - 1
                  IF( IWORK( K ).LT.0 ) THEN
                     CALL DSWAP( P, C( 1, J-IWS+1 ), 1, C( 1, K-IWS+1 ),
     $                           1 )
                     IWORK( K ) = -IWORK( K )
                     J = K
                     GO TO 60
                  END IF
               END IF
   70       CONTINUE
         END IF
      END IF
C
C                                        ( R11 R12 )
C     Logically partition E(1:NR,1:NR) = (         ),
C                                        (  0  R22 )
C
C     where R11 = R(1:RANKE,1:RANKE).
C
C     Neglect R22 and annihilate R12.
C
      IF( RANKE.LT.NR ) THEN
C
C        The block R22 is non-empty, further computations required.
C
C        Workspace: need   2*RANKE,
C                   prefer RANKE+RANKE*NB.
C
         CALL DTZRZF( RANKE, NR, E, LDE, DWORK, DWORK( RANKE+1 ),
     $                LDWORK-RANKE, INFO )
         MAXWRK = MAX( MAXWRK, INT( DWORK( RANKE+1 ) ) + RANKE )
C
C        Perform the same transformations on A(1:NR,1:NR).
C
C        Workspace: need   RANKE+NR,
C                   prefer RANKE+NR*NB.
C
         NA = NR - RANKE
         CALL DORMRZ( 'Right', 'Transpose', NR, NR, RANKE, NA, E, LDE,
     $                DWORK, A, LDA, DWORK( RANKE+1 ), LDWORK-RANKE,
     $                INFO )
         MAXWRK = MAX( MAXWRK, INT( DWORK( RANKE+1 ) ) + RANKE )
C
         IF( LRUPD ) THEN
C
C           Update C.
C
C           Workspace: need   RANKE+P,
C                      prefer RANKE+P*NB.
C
            CALL DORMRZ( 'Right', 'Transpose', P, NR, RANKE, NA, E, LDE,
     $                   DWORK, C, LDC, DWORK( RANKE+1 ), LDWORK-RANKE,
     $                   INFO )
            MAXWRK = MAX( MAXWRK, INT( DWORK( RANKE+1 ) ) + RANKE )
         END IF
C
C        Step 3: Determine the rank of A(RANKE+1:NR,RANKE+1:NR).
C
         DO 80 I = IWS, IWS + NA - 1
            IWORK( I ) = 0
   80    CONTINUE
C
C        Perform a rank-revealing QR decomposition of
C        A(RANKE+1:NR,RANKE+1:NR) with respect to the tolerance TOLDEF.
C
         SVLMAX = DLANGE( 'Frobenius', NR, NR, A, LDA, DWORK )
C
         IF( LRUPD ) THEN
C
C           IF JOBEIG = 'A' or UPDATE = 'U', then copy
C           A(RANKE+1:NR,RANKE+1:NR) to DWORK and perform the operations
C           on DWORK.
C
C           Workspace: need   NA*NA+4*NA+4,
C                      prefer NA*NA+3*NA+( NA+1 )*NB+3.
C
            ITAU = NA*NA + 1
            ISV  = ITAU  + NA
            IWRK = ISV   + 3
            CALL DLACPY( 'Full', NA, NA, A( RANKE+1, RANKE+1 ), LDA,
     $                   DWORK, NA )
            CALL MB03OD( 'QR Decomposition', NA, NA, DWORK, NA,
     $                   IWORK( IWS ), TOLDEF, SVLMAX, DWORK( ITAU ),
     $                   RANKA, DWORK( ISV ), DWORK( IWRK ),
     $                   LDWORK-IWRK+1, INFO )
         ELSE
C
C           Workspace: need   4*NA+4,
C                      prefer 3*NA+( NA+1 )*NB+3.
C
            ISV  = NA  + 1
            IWRK = ISV + 3
            CALL MB03OD( 'QR Decomposition', NA, NA,
     $                   A( RANKE+1, RANKE+1 ), LDA, IWORK( IWS ),
     $                   TOLDEF, SVLMAX, DWORK, RANKA, DWORK( ISV ),
     $                   DWORK( IWRK ), LDWORK-IWRK+1, INFO )
         END IF
C
         MAXWRK = MAX( MAXWRK, INT( DWORK( IWRK ) ) + IWRK - 1 )
C
C        Check whether there are singular values close to the
C        tolerance TOLDEF.
C
         IF( MIN( ABS( DWORK( ISV+1 ) - TOLDEF ),
     $            ABS( DWORK( ISV+2 ) - TOLDEF ) ).LT.TOLDEF/TEN )
C
C           One singular value is close to TOLDEF, warning returned.
C
     $      IWARN = 1
C
C        Set R12 to 0.
C
         IF( NA.GT.0 )
     $      CALL DLASET( 'Full', RANKE, NA, ZERO, ZERO, E( 1, RANKE+1 ),
     $                   LDE )
      ELSE
C
C        The block R22 is empty, no computation required.
C
         RANKA = 0
         NA    = 0
      END IF
C
C     Set the lower triangle of T in (1) to 0.
C
      IF( NR.GT.1 )
     $   CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, E( 2, 1 ), LDE )
C
C     Step 4: Determine if the transfer function is proper.
C
      IF( NA.EQ.RANKA ) THEN
C
C        The transfer function is proper.
C
         AB13ID = .TRUE.
      ELSE
C
C        The transfer function is improper.
C
         AB13ID = .FALSE.
      END IF
C
      DWORK( 1 ) = MAXWRK
      RETURN
C *** Last line of AB13ID ***
      END

ext: f

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
program main
  use fygg

  character(len=32) :: arg
  integer :: msg_count, msg_size
  type(yggcomm) :: outq
  logical :: ret
  integer :: i, count
  integer :: exit_code = 0
  character(len=:), allocatable :: test_msg

  if (command_argument_count().ne.2) then
     write(*, '("Error in Fortran pipe_src: The message count and size &
          &must be provided as input arguments.")')
     stop 1
  end if

  call get_command_argument(1, arg)
  read(arg, *) msg_count
  call get_command_argument(2, arg)
  read(arg, *) msg_size
  write(*, '("Hello from Fortran pipe_src: msg_count = ",i5.1,", &
       &msg_size = ",i5.1)') msg_count, msg_size

  ! Ins/outs matching with the the model yaml
  outq = ygg_output("output_pipe")
  write(*, '("pipe_src(F): Created I/O channels")')

  ! Create test message
  allocate(character(len=(msg_size + 1)) :: test_msg)
  do i = 1, msg_size
     test_msg(i:i) = '0'
  end do
  test_msg((msg_size+1):(msg_size+1)) = c_null_char

  ! Send test message multiple times
  count = 0
  do i = 1, msg_count
     ret = ygg_send(outq, test_msg, msg_size)
     if (.not.ret) then
        write(*, '("pipe_src(F): SEND ERROR ON MSG ",i5.1)') i
        exit_code = -1
        exit
     end if
     count = count + 1
  end do

  write(*, '("Goodbye from Fortran source. Sent ",i5.1," messages.")') count
  if (allocated(test_msg)) deallocate(test_msg)
  if (exit_code.lt.0) then
     stop 1
  end if

end program main

max_stars_repo_name: Xyzic/yggdrasil

max_issues_repo_name: Xyzic/yggdrasil

max_forks_repo_name: Xyzic/yggdrasil


text:
!--------------------------------------------------------------------------------------------------
! MODULE: system_var
!> @author Marin Sapunar, Ruđer Bošković Institute
!> @date October, 2017
!
! DESCRIPTION:
!> @brief Defines type for holding all variables about the system in the dynamics calculation.
!--------------------------------------------------------------------------------------------------
module system_var
    ! Import variables
    use global_defs
    use constants
    implicit none

    private
    public :: t, memory
    public :: trajtype
    public :: trajectory_next
    public :: trajectory_rewind
    public :: trajectory_write_backup
    public :: trajectory_read_backup

    !----------------------------------------------------------------------------------------------
    ! TYPE: TrajType
    !> @brief Holds dynamics variables.
    !----------------------------------------------------------------------------------------------
    type trajtype
        integer :: step = -1 !< Current step.
        integer :: substep = -1 !< Signal that substeps are currently being performed.
        real(dp) :: time = -1.0_dp !< Current time.
        real(dp) :: last_hop = -1.0_dp !< Time at which the last hop occured.

        ! Full system:
        integer :: ndim = 3 !< Number of dimensions.
        integer :: natom = 0 !< Total number of atoms.
        character(len=2), allocatable :: sym(:) !< List of atom symbols.
        real(dp), allocatable :: chrg(:) !< List of atom partial charges.
        real(dp), allocatable :: mass(:) !< List of atom masses.
        real(dp), allocatable :: geom(:, :) !< Geometry.
        real(dp), allocatable :: velo(:, :) !< Velocity.
        real(dp), allocatable :: grad(:, :) !< Gradient.

        ! QM system:
        integer :: qnatom = 0 !< Number of QM atoms.
        integer, allocatable :: qind(:) !< Indexes of the QM atoms in the full system.
        integer :: nstate = 0
        integer :: nspin = 1 !< Number of different spins. 1 or 2 supported.
        integer :: nstate_spin(2) = 0 !< Number of singlet/triplet states.
        integer :: cstate = 0 !< Current state.
      ! integer :: cspin = 1 !< Current spin.
        real(dp), allocatable :: qe(:) !< Potential energies of the QM system.
        real(dp), allocatable :: qo(:) !< Oscillator strengths of the QM system.

        ! MM system:
        integer :: mnatom = 0 !< Number of MM atoms.
        integer, allocatable :: mind(:) !< Indexes of the MM atoms in the full system.
        real(dp) :: me(2) = 0.0_dp !< MM energies (1 - QM model sys, 2 - MM sys).

        ! Surface hopping variables:
        complex(dp), allocatable :: cwf(:) !< Coeffs of the el. states in the total wf.
        real(dp), allocatable :: prob(:) !< Probabilities of hopping during current step.
        real(dp), allocatable :: olap(:, :) !< Overlaps between wfs between this and previous step.
        real(dp), allocatable :: nadv(:, :, :) !< Nonadiabatic coupling vectors.
        real(dp), allocatable :: soc(:, :) !< Spin orbit couplings.

        real(dp) :: pbcbox(6) = 0.0_dp
    contains
        procedure :: writestep => traj_writestep
        procedure :: writeheader => traj_writeheader
        procedure :: tote => traj_tote !< Total energy of the system.
        procedure :: kine => traj_kine !< Kinetic energy of the system.
        procedure :: pote => traj_pote !< Potential energy of the system.
        procedure :: qkine => traj_qkine !< QM kinetic energy.
        procedure :: mkine => traj_mkine !< MM kinetic energy.
    end type trajtype

    integer, parameter :: memory = 5
    type(trajtype), allocatable :: t(:)

contains

    !----------------------------------------------------------------------------------------------
    ! FUNCTION: ekin
    !> @brief Calculate the kinetic energy for a set of particles.
    !----------------------------------------------------------------------------------------------
    pure function ekin(mass, velo) result (en)
        real(dp), intent(in) :: mass(:)
        real(dp), intent(in) :: velo(:, :)
        real(dp) :: en
        integer :: i

        en = 0.0_dp
        do i = 1, size(mass)
            en = en + mass(i) * dot_product(velo(:, i), velo(:, i)) * 0.5_dp
        end do
    end function ekin

    pure function traj_tote(t) result(tote)
        class(trajtype), intent(in) :: t
        real(dp) :: tote
        tote = t%pote() + t%kine()
    end function traj_tote

    pure function traj_pote(t) result(pote)
        class(trajtype), intent(in) :: t
        real(dp) :: pote
        pote = t%qe(t%cstate) + t%me(2) - t%me(1)
    end function traj_pote

    pure function traj_kine(t) result(kine)
        class(trajtype), intent(in) :: t
        real(dp) :: kine
        kine = ekin(t%mass, t%velo)
    end function traj_kine

    pure function traj_qkine(t) result(qkine)
        class(trajtype), intent(in) :: t
        real(dp) :: qkine
        qkine = ekin(t%mass(t%qind), t%velo(:, t%qind))
    end function traj_qkine

    pure function traj_mkine(t) result(mkine)
        class(trajtype), intent(in) :: t
        real(dp) :: mkine
        mkine = ekin(t%mass(t%mind), t%velo(:, t%mind))
    end function traj_mkine

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: trajectory_next
    !> @brief Move array of trajectory values forward by one step at position pos.
    !----------------------------------------------------------------------------------------------
    subroutine trajectory_next(pos, t, dt, increment_time)
        integer, intent(in) :: pos
        type(trajtype), allocatable :: t(:)
        real(dp), intent(in) :: dt
        logical, intent(in) :: increment_time
        integer :: i

        do i = size(t), pos+1, -1
            t(i) = t(i-1)
        end do
        t(pos)%step = t(pos+1)%step + 1
        if (increment_time) then
            t(pos)%time = t(pos)%time + dt
        end if
    end subroutine trajectory_next

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: trajectory_rewind
    !> @brief Move array of trajectory values forward by one or more steps.
    !----------------------------------------------------------------------------------------------
    subroutine trajectory_rewind(t, nstep, increment_step)
        type(trajtype), allocatable :: t(:)
        integer, intent(in) :: nstep
        logical :: increment_step !< Increment step instead of copying value.
        integer :: i, step

        step = t(1)%step + 1
        do i = 1, size(t) - nstep
            t(i) = t(i+nstep)
        end do
        do i = size(t) - nstep + 1, size(t)
            if (allocated(t(i)%geom)) deallocate(t(i)%geom)
            if (allocated(t(i)%velo)) deallocate(t(i)%velo)
            if (allocated(t(i)%grad)) deallocate(t(i)%grad)
            if (allocated(t(i)%olap)) deallocate(t(i)%olap)
            if (allocated(t(i)%nadv)) deallocate(t(i)%nadv)
        end do
        if (increment_step) t(1)%step = step
    end subroutine trajectory_rewind

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: Traj_WriteStep
    !
    ! DESCRIPTION:
    !> @brief Write requested output from each step.
    !----------------------------------------------------------------------------------------------
    subroutine traj_writestep(t, popt, res_dir)
        class(trajtype) :: t
        logical, intent(in) :: popt(:) !< Print options.
        character(len=*), intent(in) :: res_dir ! Directory for writing results.
        integer :: ounit, i
        real(dp) :: time_fs

        time_fs = t%time * aut_fs
        if (popt(1)) then
            open(newunit=ounit, file=res_dir//'/energy.dat', action='write', position='append')
            write(ounit, 1001, advance='no') time_fs, t%cstate
            write(ounit, 1002, advance='no') t%tote(), t%qe(t%cstate)
            write(ounit, 1002, advance='no') t%qe(:)
            write(ounit, *)
            close(ounit)

            if (t%mnatom > 0) then
                open(newunit=ounit, file=res_dir//'/mm.dat', action='write', position='append')
                write(ounit, 1001, advance='no') time_fs
                write(ounit, 1002, advance='no') t%mkine(), t%me
                write(ounit, 1002, advance='no') t%pbcbox
                write(ounit, *)
                close(ounit)
            end if
        end if

        if (popt(2)) then
            open(newunit=ounit, file=res_dir//'/trajectory.xyz', action='write', position='append')
            write(ounit, *) t%natom
            write(ounit, *) 't= ', time_fs, 'fs, state=', t%cstate
            do i = 1, t%natom
                write(ounit, 1003) t%sym(i), t%geom(:, i) * a0_A
            end do
            close(ounit)
        end if

        if (popt(3)) then
            open(newunit=ounit, file=res_dir//'/geometry', action='write', position='append')
            do i = 1, t%natom
                write(ounit, 1002) t%geom(:, i)
            end do
            close(ounit)
        end if

        if (popt(4)) then
            open(newunit=ounit, file=res_dir//'/velocity', action='write', position='append')
            do i = 1, t%natom
                write(ounit, 1002) t%velo(:, i)
            end do
            close(ounit)
        end if

        if (popt(5)) then
            open(newunit=ounit, file=res_dir//'/gradient', action='write', position='append')
            do i = 1, t%natom
                write(ounit, 1002) t%grad(:, i)
            end do
            close(ounit)
        end if

        if (popt(6)) then
            open(newunit=ounit, file=res_dir//'/cwf.dat', action='write', position='append')
            write(ounit, 1006) abs(t%cwf)
            close(ounit)
        end if

      ! if (popt(8)) then
      !     open(newunit=ounit, file=res_dir//'/prob.dat', action='write', position='append')
      !     write(ounit, 1006) t%prob
      !     close(ounit)
      ! end if

        if (popt(7)) then
            open(newunit=ounit, file=res_dir//'/overlap', action='write', position='append')
            write(ounit, *) 't= ', time_fs, 'fs, state=', t%cstate
            do i = 1, t%nstate
                write(ounit, 1006) t%olap(i, 1:t%nstate)
            end do
            close(ounit)
        end if

        if (popt(10)) then
            open(newunit=ounit, file=res_dir//'/oscill.dat', action='write', position='append')
            write(ounit, 1006) t%qo
            close(ounit)
        end if

        if (popt(11)) then
            open(newunit=ounit, file=res_dir//'/qm_traj.xyz', action='write', position='append')
            write(ounit, *) t%qnatom
            write(ounit, *) 't= ', time_fs, 'fs, state=', t%cstate
            do i = 1, t%qnatom
                write(ounit, 1003) t%sym(t%qind(i)), t%geom(:, t%qind(i)) * a0_A
            end do
            close(ounit)
        end if

1001 format (f12.5,2x,i3)
1002 format (1x,1000f18.10)
1003 format (1x,a2,2x,1000f18.10)
1006 format (1x,1000e22.12)
    end subroutine traj_writestep

    subroutine trajectory_write_backup(backup_file, t)
        character(len=*), intent(in) :: backup_file
        class(trajtype), intent(in) :: t(:)
        integer :: i, ounit

        open(newunit=ounit, file=backup_file, action='write')
        write(ounit, *) memory
        do i = 1, memory
            write(ounit, *) t(i)%step
            write(ounit, *) t(i)%substep
            write(ounit, *) t(i)%time
            write(ounit, *) t(i)%ndim
            write(ounit, *) t(i)%natom
            write(ounit, *) allocated(t(i)%sym)
            if (allocated(t(i)%sym)) write(ounit, *) t(i)%sym(:)
            write(ounit, *) allocated(t(i)%chrg)
            if (allocated(t(i)%chrg)) write(ounit, *) t(i)%chrg(:)
            write(ounit, *) allocated(t(i)%mass)
            if (allocated(t(i)%mass)) write(ounit, *) t(i)%mass(:)
            write(ounit, *) allocated(t(i)%geom)
            if (allocated(t(i)%geom)) write(ounit, *) t(i)%geom(:, :)
            write(ounit, *) allocated(t(i)%velo)
            if (allocated(t(i)%velo)) write(ounit, *) t(i)%velo(:, :)
            write(ounit, *) allocated(t(i)%grad)
            if (allocated(t(i)%grad)) write(ounit, *) t(i)%grad(:, :)
            write(ounit, *) t(i)%qnatom
            write(ounit, *) allocated(t(i)%qind)
            if (allocated(t(i)%qind)) write(ounit, *) t(i)%qind(:)
            write(ounit, *) t(i)%nspin
            write(ounit, *) t(i)%nstate
            write(ounit, *) t(i)%nstate_spin
            write(ounit, *) t(i)%cstate
            write(ounit, *) allocated(t(i)%qe)
            if (allocated(t(i)%qe)) write(ounit, *) t(i)%qe(:)
            write(ounit, *) allocated(t(i)%qo)
            if (allocated(t(i)%qo)) write(ounit, *) t(i)%qo(:)
            write(ounit, *) t(i)%mnatom
            write(ounit, *) allocated(t(i)%mind)
            if (allocated(t(i)%mind)) write(ounit, *) t(i)%mind(:)
            write(ounit, *) t(i)%me(2)
            write(ounit, *) allocated(t(i)%cwf)
            if (allocated(t(i)%cwf)) write(ounit, *) t(i)%cwf(:)
            write(ounit, *) allocated(t(i)%prob)
            if (allocated(t(i)%prob)) write(ounit, *) t(i)%prob(:)
            write(ounit, *) allocated(t(i)%olap)
            if (allocated(t(i)%olap)) write(ounit, *) t(i)%olap(:, :)
            write(ounit, *) allocated(t(i)%nadv)
            if (allocated(t(i)%nadv)) write(ounit, *) t(i)%nadv(:, :, :)
            write(ounit, *) t(i)%pbcbox(6)
        end do
        close(ounit)
    end subroutine trajectory_write_backup

    subroutine trajectory_read_backup(backup_file, t)
        character(len=*), intent(in) :: backup_file
        class(trajtype), intent(inout) :: t(:)
        integer :: i, iunit
        logical :: check

        open(newunit=iunit, file=backup_file, status='old', action='read')
        read(iunit, *) i
        if (i /= memory) stop 'Error. Bad memory parameter in system_var.'
        do i = 1, memory
            read(iunit, *) t(i)%step
            read(iunit, *) t(i)%substep
            read(iunit, *) t(i)%time
            read(iunit, *) t(i)%ndim
            read(iunit, *) t(i)%natom
            read(iunit, *) check
            if (check) then
                allocate(t(i)%sym(t(i)%natom))
                read(iunit, *) t(i)%sym(:)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%chrg(t(i)%natom))
                read(iunit, *) t(i)%chrg(:)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%mass(t(i)%natom))
                read(iunit, *) t(i)%mass(:)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%geom(t(i)%ndim, t(i)%natom))
                read(iunit, *) t(i)%geom(:, :)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%velo(t(i)%ndim, t(i)%natom))
                read(iunit, *) t(i)%velo(:, :)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%grad(t(i)%ndim, t(i)%natom))
                read(iunit, *) t(i)%grad(:, :)
            end if
            read(iunit, *) t(i)%qnatom
            read(iunit, *) check
            if (check) then
                allocate(t(i)%qind(t(i)%qnatom))
                read(iunit, *) t(i)%qind(:)
            end if
            read(iunit, *) t(i)%nspin
            read(iunit, *) t(i)%nstate
            read(iunit, *) t(i)%nstate_spin
            read(iunit, *) t(i)%cstate
            read(iunit, *) check
            if (check) then
                allocate(t(i)%qe(t(i)%nstate))
                read(iunit, *) t(i)%qe(:)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%qo(t(i)%nstate_spin(1)-1))
                read(iunit, *) t(i)%qo(:)
            end if
            read(iunit, *) t(i)%mnatom
            read(iunit, *) check
            if (check) then
                allocate(t(i)%mind(t(i)%mnatom))
                read(iunit, *) t(i)%mind(:)
            end if
            read(iunit, *) t(i)%me(2)
            read(iunit, *) check
            if (check) then
                allocate(t(i)%cwf(t(i)%nstate))
                read(iunit, *) t(i)%cwf(:)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%prob(t(i)%nstate))
                read(iunit, *) t(i)%prob(:)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%olap(t(i)%nstate, t(i)%nstate))
                read(iunit, *) t(i)%olap(:,:)
            end if
            read(iunit, *) check
            if (check) then
                allocate(t(i)%nadv(t(i)%ndim*t(i)%natom, t(i)%nstate, t(i)%nstate))
                read(iunit, *) t(i)%nadv(:, :, :)
            end if
            read(iunit, *) t(i)%pbcbox(6)
        end do
        close(iunit)
    end subroutine trajectory_read_backup

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: Traj_writeheader
    !
    ! DESCRIPTION:
    !> @brief Initialize output files and write column headers.
    !----------------------------------------------------------------------------------------------
    subroutine traj_writeheader(t, popt, res_dir)
        class(trajtype) :: t
        logical, intent(in) :: popt(:) !< Print options.
        character(len=*), intent(in) :: res_dir ! Directory for writing results.
        integer :: ounit, i
        real(dp) :: time_fs

        time_fs = t%time * aut_fs
        if (popt(1)) then
            open(newunit=ounit, file=res_dir//'/energy.dat', action='write', position='append')
            write(ounit, '(a)', advance='no') '#'
            write(ounit, '(a12,1x,a4)', advance='no') 't,', 'cst,'
            write(ounit, '(1x,2a18)', advance='no') 'etot,', 'epot,'
            write(ounit, '(1x)', advance='no')
            do i = 1, size(t%qe)
                write(ounit, '(12x,a2,i3.3)', advance='no') 'qe', i
                if (i < size(t%qe)) write(ounit, '(a)', advance='no') ','
            end do
            write(ounit, *)
            close(ounit)

            if (t%mnatom > 0) then
                open(newunit=ounit, file=res_dir//'/mm.dat', action='write', position='append')
                write(ounit, '(a)', advance='no') '#'
                write(ounit, '(a12,1x)', advance='no') 't,'
                write(ounit, '(2x,3a18)', advance='no') 'mekin,', 'me1,', 'me2,'
                write(ounit, '(1x,5a18,a17)', advance='no') 'pbc1,', 'pbc2,', 'pbc3,', 'pbc4,', 'pbc5,', 'pbc6'
                write(ounit, *)
                close(ounit)
            end if
        end if
    end subroutine traj_writeheader

end module system_var

max_stars_repo_name: marin-sapunar/ZagHop

max_issues_repo_name: marin-sapunar/ZagHop

max_forks_repo_name: marin-sapunar/ZagHop