c Copyright (C) 2008-2014 Vincent Eymet c c KDISTRIBUTION is free software; you can redistribute it and/or modify c it under the terms of the GNU General Public License as published by c the Free Software Foundation; either version 3, or (at your option) c any later version. c KDISTRIBUTION is distributed in the hope that it will be useful, c but WITHOUT ANY WARRANTY; without even the implied warranty of c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the c GNU General Public License for more details. c You should have received a copy of the GNU General Public License c along with KDISTRIBUTION; if not, see c program main implicit none include 'max.inc' include 'formats.inc' include 'mpif.h' integer kind,kpts,i,j,Nb,band,imin,quad,Nq,m,ns integer Np,nk,i_min,i_max,nl,valid,ik,mol integer band1,band2 double precision alt(1:Nmax) double precision pres(1:Nmax) double precision temp(1:Nmax) double precision x(1:Nmax,1:Nmol_mx) character*(Nchar_mx) filelist(1:Nmax) double precision eps1,eps2 double precision alpha,beta,b(Nqmx) double precision endpts(2),fixed_g(2) double precision g(Nqmx) double precision w(Nqmx) double precision sum_w double precision nu_min(1:Nbmx),nu_max(1:Nbmx) double precision nu(1:kmax),k(1:kmax) double precision z(1:kmax) double precision k_min,k_max double precision gac(0:Npmax),kac(0:Npmax) double precision epsilon,raf integer sit double precision err(1:Nmax,1:Nbmx,1:Nqmx) double precision kdist(1:Nmax,1:Nbmx,1:Nqmx) double precision delta_k(1:Nmax,1:Nbmx,1:Nqmx) double precision tr_exact(1:Nmax,1:Nbmx) double precision tr_quad(1:Nmax,1:Nbmx) double precision tr_exactp(1:Nbmx) double precision tr_quadp(1:Nbmx) double precision tr_exact0(1:Nbmx) double precision tr_quad0(1:Nbmx) double precision kdistp(1:Nbmx,1:Nqmx) double precision delta_kp(1:Nbmx,1:Nqmx) double precision errp(1:Nbmx,1:Nqmx) double precision kdist0(1:Nbmx,1:Nqmx) double precision delta_k0(1:Nbmx,1:Nqmx) double precision err0(1:Nbmx,1:Nqmx) double precision g_step(1:Nsmx) character*(Nchar_mx) kfile,base,quad_file logical debug logical k_is_null integer ios,iosP,iosT,iosX integer Nlines,Nvalues,Nvalues_read integer Nlines_to_first_nu,Nlines_preamble double precision nu_prev,k_prev double precision nu_tmp,k_tmp integer Nmol(1:Nmax),Nmolec character*(Nchar_mx) molec_names(1:Nmax,1:Nmol_mx) integer ilev,imol c MPI integer ierr,errorcode,nproc,pindex integer stat(MPI_STATUS_SIZE) integer chunk,nchunks integer code0,codep integer proc,fin,fpf,nsp integer free(1:Nproc_mx) integer ch_idx(1:Nproc_mx) c MPI c label integer strlen character*(Nchar_mx) label label='program main' c MPICH call MPI_INIT(ierr) if (ierr.NE.MPI_SUCCESS) THEN write(*,*) 'Error from program ', & label(1:strlen(label)),' :' write(*,*) 'MPI could not be initialized' call MPI_ABORT(MPI_COMM_WORLD,errorcode,ierr) endif c Get 'np', the number of processes call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ierr) if (nproc.lt.2) then write(*,*) 'Error from program ', & label(1:strlen(label)),' :' write(*,*) 'Number of processes is:',nproc write(*,*) 'and it should be greater than 1' stop endif c Get 'pindex', the index of the current process call MPI_COMM_RANK(MPI_COMM_WORLD,pindex,ierr) c MPICH if (pindex.eq.0) then ! root process only write(*,*) 'KDISTRIBUTION Copyright (C) 2008-2014 V. Eymet' write(*,*) 'This program comes with ABSOLUTELY NO WARRANTY;' write(*,*) 'This is free software, and you are welcome to'// & ' redistribute it' write(*,*) 'under certain conditions; see file "COPYING"'// & ' for details.' write(*,*) debug=.false. c ********** FIRST STEP ************** c production of quadrature abscissas and weights c base='./data/hires_spectra' write(*,*) 'Production of quadrature abscissas and weights...' c Reading options c warning: only Legendre quadrature (kind=1 and kpts=0), c Radau (kind=1 and kpts=1) and Lobatto (kind=1 and kpts=2) are available c Debug if (debug) then write(*,*) 'This is: ',label(1:strlen(label)) write(*,*) 'Preparing to read options.in' endif c Debug call read_options(kind,kpts,Nq,fixed_g,alpha,beta, & endpts,Np,epsilon,sit,quad_file) c Debug if (debug) then write(*,*) 'This is: ',label(1:strlen(label)) write(*,*) 'options.in has been read' write(*,*) 'Now preparing to compute quadrature abscissas' write(*,*) 'with input parameter kind=',kind endif c Debug if (kind.eq.0) then c reading quadrature abscissas g(quad) and weights w(quad) open(10,file=quad_file(1:strlen(quad_file))) read(10,*) do quad=1,Nq read(10,*) g(quad),w(quad) enddo close(10) c Debug c do quad=1,Nq c write(*,*) 'g(',quad,')=',g(quad),w(quad) c enddo c Debug else c producing quadrature abscissas t(quad) and weights w1(quad) if (kind.eq.1) then ! Legendre call myquad(Nq,1,alpha,beta,g,w,valid) else if (kind.eq.2) then ! Legendre-Radau call myquad(Nq,2,alpha,beta,g,w,valid) else if (kind.eq.3) then ! Legendre-Lobatto call myquad(Nq,3,alpha,beta,g,w,valid) else if (kind.eq.4) then ! Chebyshev first kind call gaussq(2,Nq,alpha,beta,kpts,endpts,b,g,w) else if (kind.eq.5) then ! Chebyshev second kind call gaussq(3,Nq,alpha,beta,kpts,endpts,b,g,w) else if (kind.eq.6) then ! Hermite call gaussq(4,Nq,alpha,beta,kpts,endpts,b,g,w) else if (kind.eq.7) then ! Jacobi call myquad(Nq,4,alpha,beta,g,w,valid) else if (kind.eq.8) then ! Laguerre call gaussq(6,Nq,alpha,beta,kpts,endpts,b,g,w) endif endif c Debug if (debug) then write(*,*) 'This is: ',label(1:strlen(label)) write(*,*) 'Quadrature abscissas have been computed' endif c Debug c sort abscissas and weights by ascending order call sort_gw(Nq,g,w) c Debug if (debug) then write(*,*) 'This is: ',label(1:strlen(label)) write(*,*) 'Quadrature abscissas and weights ' & //'have been sorted' endif c Debug c writing output information call write_quad_results(kpts,kind,Nq,fixed_g,g,w) c Debugâ—Š if (debug) then write(*,*) 'This is: ',label(1:strlen(label)) write(*,*) 'Quadrature results have been recorded' endif c Debug write(*,*) '...done' c c ************************************ c ********** SECOND STEP ************** c Production of the k-distributions write(*,*) 'Production of k-distribution data set...' call make_model(base,m,alt,pres,temp, & Nmol,molec_names,x,filelist) c Debug c write(*,*) 'In main: m=',m c stop c do ilev=1,m c write(*,*) 'level:',ilev c write(*,*) 'alt=',alt(ilev) c write(*,*) 'pres=',pres(ilev) c write(*,*) 'temp=',temp(ilev) c write(*,*) 'Nmol=',Nmol(ilev) c write(*,*) 'x=',(x(ilev,imol),imol=1,Nmol(ilev)) c enddo ! ilev c stop c Debug c looking for values of k(g(i)) call read_nu(Nb,nu_min,nu_max) call read_data(i_min,i_max) call look_for_errors(m,i_min,i_max) nchunks=i_max-i_min+1 code0=1 ! data transfer do proc=1,nproc-1 c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'sending code0=',code0,' to proc=',proc endif c Debug call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'code0 has been sent' endif c Debug enddo c send integer data c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'sending data through MPI_BCAST' endif c Debug call MPI_BCAST(debug,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nb,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nq,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Np,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nmol,Nmax,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sit,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) c send double precision data call MPI_BCAST(epsilon,1, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(eps1,1, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(eps2,1, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nu_min,Nbmx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nu_max,Nbmx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(w,Nqmx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(g,Nqmx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(alt,Nmax, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(pres,Nmax, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(temp,Nmax, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(x,Nmax*Nmol_mx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(filelist,Nmax*Nchar_mx, & MPI_CHAR & ,0,MPI_COMM_WORLD,ierr) c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'all data has been sent' endif c Debug do proc=1,nproc-1 free(proc)=1 enddo nk=0 fin=0 chunk=0 do while (fin.eq.0) 444 continue call find_free_process(free,nproc-1,fpf,proc) if (fpf.eq.1) then ! a free child process found chunk=chunk+1 if (chunk.le.nchunks) then ch_idx(proc)=chunk code0=2 c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'sending code0=',code0,' to proc=',proc endif c Debug call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) i=i_min+chunk-1 c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'sending i=',i,' to proc=',proc endif c Debug call MPI_SEND(i,1,MPI_INTEGER, & proc,1,MPI_COMM_WORLD,ierr) free(proc)=0 ! child process index 'proc' is busy else ! chunk>nchunks free(proc)=2 ! means child process index 'proc' has been stopped call MPI_SEND(0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) endif goto 444 else ! no free child process found call stopped_processes(free,nproc-1,nsp) if (nsp.eq.nproc-1) then fin=1 goto 111 endif c wait for results c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'waiting data' endif c Debug call MPI_RECV(proc,1,MPI_INTEGER & ,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD, & stat,ierr) c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'proc=',proc endif c Debug call MPI_RECV(kdist0,Nbmx*Nqmx,MPI_DOUBLE_PRECISION & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) call MPI_RECV(delta_k0,Nbmx*Nqmx,MPI_DOUBLE_PRECISION & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) call MPI_RECV(err0,Nbmx*Nqmx,MPI_DOUBLE_PRECISION & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) call MPI_RECV(tr_exact0,Nbmx,MPI_DOUBLE_PRECISION & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) call MPI_RECV(tr_quad0,Nbmx,MPI_DOUBLE_PRECISION & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) free(proc)=1 ! child process index 'proc' is free again i=i_min+ch_idx(proc)-1 ! index of current level c Debug c do quad=1,Nq c write(*,*) quad,kdist0(1,quad) c enddo c Debug c record results for the current level call write_kdist_results_tmp(i,Nb,Nq,nu_min,nu_max, & w,kdist0,delta_k0,tr_exact0,tr_quad0, & pres,temp,Nmol,molec_names,x, & filelist(i)(1:strlen(filelist(i)))) do band=1,Nb do quad=1,Nq kdist(i,band,quad)=kdist0(band,quad) delta_k(i,band,quad)=delta_k0(band,quad) err(i,band,quad)=err0(band,quad) enddo ! quad tr_exact(i,band)=tr_exact0(band) tr_quad(i,band)=tr_quad0(band) enddo ! band endif ! fpf=0 or 1 111 continue enddo ! while fin=0 c$$$c stopping all child processes c$$$ do proc=1,nproc-1 c$$$ code0=0 c$$$c Debug c$$$ if (debug) then c$$$ write(*,*) 'This is process:',pindex c$$$ write(*,*) 'sending code0=',code0,' to proc=',proc c$$$ endif c$$$c Debug c$$$ call MPI_SEND(code0,1,MPI_INTEGER c$$$ & ,proc,1,MPI_COMM_WORLD,ierr) c$$$ enddo ! proc c ************************************ c writing output information c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'writing kdist.txt' endif c Debug call write_kdist_results(i_min,i_max,Nb,Nq,nu_min,nu_max,w, & kdist,delta_k,tr_exact,tr_quad, & pres,temp,Nmol,molec_names,x,filelist) write(*,*) '...done' c ************************************ write(*,*) 'cleaning temporary files...' call clean() write(*,*) '...done' else ! child processes 333 continue c wait for "codep" from root process c Debug c if (debug) then c write(*,*) 'This is process:',pindex c write(*,*) 'waiting for codep' c endif c Debug call MPI_RECV(codep,1,MPI_INTEGER & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) c Debug c if (debug) then c write(*,*) 'This is process:',pindex c write(*,*) 'data received: codep=',codep c endif c Debug c codep=0 means exit if (codep.eq.0) then goto 666 c codep=1 means reception of data else if (codep.eq.1) then c Debug c if (debug) then c write(*,*) 'This is process:',pindex c write(*,*) 'waiting for data' c endif c Debug c send integer data call MPI_BCAST(debug,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nb,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nq,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Np,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nmol,Nmax,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sit,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) c send double precision data call MPI_BCAST(epsilon,1, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(eps1,1, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(eps2,1, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nu_min,Nbmx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nu_max,Nbmx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(w,Nqmx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(g,Nqmx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(alt,Nmax, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(pres,Nmax, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(temp,Nmax, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(x,Nmax*Nmol_mx, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(filelist,Nmax*Nchar_mx, & MPI_CHAR & ,0,MPI_COMM_WORLD,ierr) c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'data received through MPI_BCAST' endif c Debug c codep=2 means a computation has to be performed else if (codep.eq.2) then c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'waiting for i' endif c Debug call MPI_RECV(i,1,MPI_INTEGER & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'data received: i=',i endif c Debug c Debug write(*,*) 'Computation for spectrum: ',i c Debug kfile=filelist(i)(1:strlen(filelist(i))) c opening high-resolution spectra files open(13,file=kfile(1:strlen(kfile)) & ,status='old',iostat=ios) if (ios.ne.0) then ! file not found call error(label) write(*,*) 'File could not be opened:' write(*,*) kfile(1:strlen(kfile)) stop else call get_nlines(kfile,pindex,Nlines) endif band1=1 band2=Nb do band=band1,band2 c Debug c write(*,*) 'i=',i,' band=',band c Debug if (band.eq.band1) then do j=1,5 read(13,*) enddo ! j read(13,*) Nmolec Nlines_to_first_nu=5*Nmolec+3 Nlines_preamble=Nlines_to_first_nu+6 Nvalues=Nlines-Nlines_preamble Nvalues_read=0 do j=1,Nlines_to_first_nu read(13,*) enddo ! j nu_prev=0.0D+0 k_prev=0.0D+0 nu_tmp=0.0D+0 k_tmp=0.0D+0 endif c Debug c write(*,*) 'calling read_single_spectrum' c write(*,*) 'nu_min(',band,')=',nu_min(band) c write(*,*) 'nu_max(',band,')=',nu_max(band) c Debug call read_single_kspectrum(i,band, & Nlines,Nvalues,Nvalues_read, & nu_min(band),nu_max(band),pindex, & nu_prev,k_prev, & nu_tmp,k_tmp, & nk,nu,k,Nmolec,k_is_null) c Debug c write(*,*) 'read_single_spectrum: ok' c write(*,*) 'nu(1)=',nu(1),' k(1)=',k(1) c write(*,*) 'nu(',nk,')=',nu(nk),' k(',nk,')=',k(nk) c write(*,*) 'nu(3)=',nu(3),' k(3)=',k(3) c write(*,*) 'band=',band c Debug c Debug c if (band.eq.2) then c write(*,*) 'stop required' c stop c endif c Debug if (nk.gt.kmax) then write(*,*) 'Error from program ', & label(1:strlen(label)),' :' write(*,*) 'nk=',nk,' while kmax=',kmax write(*,*) 'interface:',i,' band:',band stop endif c Debug c write(*,*) 'spectrum was read for band=',band c write(*,*) 'k_is_null=',k_is_null c Debug if (k_is_null) then do quad=1,Nq kdistp(band,quad)=0.0D+0 delta_kp(band,quad)=0.0D+0 errp(band,quad)=0.0D+0 enddo ! quad goto 412 endif if (sit.ne.1) then c Debug c write(*,*) 'ok0' c Debug call cubic_spline(nk,nu,k,z,valid) c Debug c write(*,*) 'ok1' c Debug endif call minmax(nk,k,k_min,k_max) c Debug c write(*,*) 'Now calling analyse_g' c write(*,*) 'k_min=',k_min c write(*,*) 'k_max=',k_max c write(*,*) 'valid=',valid c Debug call analyse_g(band,nk,nu,k,eps1,eps2,k_min,k_max,sit,z, & ns,g_step) c Debug c write(*,*) 'analyse_g -> ok' c Debug c production of acceleration arrays "kac" and "gac" c Debug c write(*,*) 'Now calling acceleration' c Debug call acceleration(i,band,nk,nu,k,eps1,eps2, & Np,k_max,sit,z, & kac,gac) c Debug c write(*,*) 'acceleration -> ok' c Debug c Debug c do ik=0,Np-1 c write(*,*) gac(ik),gacT(ik) c enddo c stop c Debug c Inversion of g(i) values do quad=1,Nq c Debug c write(*,*) 'calling g_inversion for quad=',quad c Debug call g_inversion(i,band,quad,g(quad),Np,gac,kac, & epsilon,nk,nu,k,eps1,eps2,sit,z, & kdistp(band,quad), & delta_kp(band,quad),errp(band,quad)) c Debug c write(*,*) 'g=',g(quad), c & ' k=',kdistp(band,quad) c write(*,*) 'kdistp(',band,',',quad,')=', c & kdistp(band,quad) c write(*,*) 'ok0' c Debug enddo ! quad c Debug c write(*,*) 'band=',band,kdistp(band,1),kdistpx(band,1,1) c Debug 412 continue c Debug c write(*,*) 'Now calling test_procedure' c Debug call test_procedure(nk,Nb,Nq,i,band,nu,k,kdistp,w, & tr_exactp(band),tr_quadp(band)) c Debug c write(*,*) 'all done for band=',band c Debug enddo ! band if (ios.eq.0) then close(13) endif c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'sending pindex=',pindex endif c Debug call MPI_SEND(pindex,1,MPI_INTEGER & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(kdistp,Nbmx*Nqmx,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(delta_kp,Nbmx*Nqmx,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(errp,Nbmx*Nqmx,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(tr_exactp,Nbmx,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(tr_quadp,Nbmx,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) c Debug if (debug) then write(*,*) 'This is process:',pindex write(*,*) 'all data was sent' endif c Debug endif ! codep c wait for a new execution code goto 333 666 continue endif ! pindex c MPICH call MPI_FINALIZE(ierr) c MPICH end