c kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum c Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet c c This file must be used under the terms of the CeCILL license. c This source file is licensed as described in the file COPYING, which c you should have received as part of this distribution. The terms c are also available at c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt c program kspectrum implicit none c ------------------------------------------------------------------------ c First, install mpich2: c retrieve the archive, unzip and untar it, install it on the system: c > ./configure --prefix=/home/[user]/parallele/mpich2 c > make c > make install c c Before running MPD: c c 1) create a ~/.mpd.conf file: c > echo secretword=[secretword] >> ~/.mpd.conf c > chmod 600 .mpd.conf c c 2) allow your machine to connect by ssh to every other machine, without password: c Let's say you want to connect from machine A to machine B without typing any password. c First, create a dsa key on machine A: c > ssh-keygen -t dsa c then leave everything blank (only use 'enter' to proceed) c Next, copy the "id_dsa.pub" file that has been generated on machine B; then c on machine B: c > cd .ssh c > cat id_dsa.pub >> authorized_keys c c 3) create the mpd.hosts file: each line must contain the name of a machine c that is part of your cluster, by order of availability: c [host1].[domain] c [host2].[domain] c [host3].[domain] c etc c ------------------------------------------------------------------------ c Running MPD: c > mpdboot -n [#] c with # the number of machines MPD has to run on c c use 'mpif77' in order to compile fortran 77 code in your makefile (or whatever) c c Code execution: c > mpirun -np [#] [executable] c with # the number of processes; in a master-slave algorithm, there is a master process and c 'np-1' slave processes. If the master process is idle most of the time, it might be a good idea c to have 'np-1' equal to the number of physical cores, so that each (time consuming) process c uses a core. In this case, 'np' is equal to the number of cores plus one. The total computation c time should be roughly divided by the number of cores (if communication times are negligeable, c and of course if the master process does not use too much time). c c In the particular case of the present code, there is no problem using np=number of cores+1 c c Exit MPD: c > mpdallexit c ------------------------------------------------------------------------ include 'mpif.h' include 'max.inc' include 'formats.inc' include 'arrays1.inc' integer l,i,j,m,molec,isotop,isof,line integer nfiles,band_min,band_max,i_min,i_max,Npass,nl_lp integer sdis,htv,htsd_co2,htsd_h2o,dr,nis logical slwr,slam,slas logical custom_db,ltr,usl,ucia,pct,rid,tfl,rl integer line_profile,sl_choice,ciac,est logical ciawr,ciatr integer uaal,uanb,sda,bslal integer dorg logical sens_T,sens_P,sens_x integer Nb,band,ik double precision sigma,k,kCIA,dkCIA_dP,dkCIA_dT,dkCIA_dxco2 double precision dsigma_dP,dsigma_dT,dsigma_dx double precision dk_dP,dk_dT,dk_dx double precision band_sigma(1:Nkmx) double precision band_k(1:Nkmx) double precision band_dsigma_dP(1:Nkmx) double precision band_dsigma_dT(1:Nkmx) double precision band_dsigma_dx(1:Nkmx,1:Nmol_max) double precision band_dk_dP(1:Nkmx) double precision band_dk_dT(1:Nkmx) double precision band_dk_dx(1:Nkmx,1:Nmol_max) double precision P(1:Nmax) double precision T(1:Nmax) double precision x(1:Nmax,1:Nmol_max) double precision Ps(1:Nmol_max) double precision nu_low(1:Nbmx) double precision nu_hi(1:Nbmx) integer nk,nsint,Nk_total,Nk_finished,Nk_remain double precision nu(1:Nkmx),mtpk double precision k_mxep,sd_mxep integer c(1:Nline_mx) logical large_db c double precision sigma_const(1:Nline_mx) c double precision k_const(1:Nline_mx) double precision dsigmac_dP double precision dsigmac_dT double precision dsigmac_dx(1:Nmol_max) double precision dkc_dP double precision dkc_dT double precision dkc_dx(1:Nmol_max) c init double precision c2,Rgp c isotopic data integer Hnmol,Nmol,Nmol_ini integer Hmol_index(1:Nmol_max) integer Hmol_niso(1:Nmol_max) character*10 Hmol_name(1:Nmol_max) integer Hiso_index(1:Nmol_max,1:Niso_max) double precision Hiso_abundance(1:Nmol_max,1:Niso_max) double precision Hiso_Qref(1:Nmol_max,1:Niso_max) integer Hiso_g(1:Nmol_max,1:Niso_max) double precision Hiso_mass(1:Nmol_max,1:Niso_max) integer mol_index(1:Nmol_max) integer mol_niso(1:Nmol_max) character*10 mol_name(1:Nmol_max) integer iso_index(1:Nmol_max,1:Niso_max) double precision abundances(1:Nmax,1:Nmol_max,1:Niso_max) double precision iso_Qref(1:Nmol_max,1:Niso_max) integer iso_g(1:Nmol_max,1:Niso_max) double precision iso_mass(1:Nmol_max,1:Niso_max) double precision iso_abundance(1:Nmol_max,1:Niso_max) integer dbr,lastf,lastl,pass,raz,lastfpp,lastlpp integer mol_index_ini(1:Nmol_max) integer mol_niso_ini(1:Nmol_max) character*10 mol_name_ini(1:Nmol_max) integer iso_index_ini(1:Nmol_max,1:Niso_max) double precision iso_Qref_ini(1:Nmol_max,1:Niso_max) integer iso_g_ini(1:Nmol_max,1:Niso_max) double precision iso_mass_ini(1:Nmol_max,1:Niso_max) c inputs of read_line_parameters character*(Nchar_mx) kfile,kfile_dP,kfile_dT,kfile_dx character*(Nchar_mx) filename(1:nfiles_mx) character*(Nchar_mx) input_file character*(Nchar_mx) planet_descriptor integer fformat(1:nfiles_mx) c outputs of read_line_parameters integer Nlines,Nlc,Nlnc integer ncl_indexes(1:Nline_mx) double precision Tref(1:Nmol_max) double precision nu_min,nu_max double precision B_mxep,constant_dnu integer nchunks double precision Sref_frac integer mol,co2f double precision T_hitemp integer i_inf integer i_sup integer band_inf integer band_sup double precision tr_dist,pr_kc,pr_sd integer ngamma integer nlw_cz,np_lc,np_bl double precision alt(1:Nmax) double precision sigmac,kc double precision numin,numax integer nk_tmp integer Nit_max,ndiv double precision eps_max,dx0,dx0_ori,pp,dp_max,y integer n1,n2 double precision xres(1:Nppc_max,1:3) c optimization integer sorted_index(1:Nline_mx) double precision nuc(1:Nline_mx) double precision dnuc_dP(1:Nline_mx) double precision gamma_L(1:Nline_mx) double precision dgammaL_dP(1:Nline_mx) double precision dgammaL_dT(1:Nline_mx) double precision dgammaL_dx(1:Nline_mx) double precision gamma_D(1:Nline_mx) double precision dgammaD_dT(1:Nline_mx) double precision S(1:Nline_mx) double precision dS_dT(1:Nline_mx) double precision densities(1:Nline_mx) double precision drho_dP(1:Nline_mx) double precision drho_dT(1:Nline_mx) double precision drho_dx(1:Nline_mx) c optimization c acceleration integer nacc,nxtab,ifile,out1,out2 double precision xtab(1:Nxtab_mx,1:3) double precision xacc(1:Nacc_mx,1:2) integer nk_func,nx(1:Nytab_mx) integer kindex(1:Nytab_mx,1:2) double precision ktab(1:Nktab_mx,1:3) double precision yvalues(1:Nytab_mx) c acceleration integer over,last_i,last_band double precision slimits1(1:Nti_mx,2) integer ncalc0,ncalc2 c MPICH integer ierr,errorcode,np,pindex,proc integer stat(MPI_STATUS_SIZE) integer Nk_p integer Nkrun(1:Nchunk_mx) integer ik1(1:Nchunk_mx) integer ik2(1:Nchunk_mx) integer free(1:Nproc_mx) integer ch_idx(1:Nproc_mx) double precision nu_tmp(1:Nkmx) double precision res(1:Nkmx,1:2) double precision dres_dP(1:Nkmx,1:2) double precision dres_dT(1:Nkmx,1:2) double precision dres_dx(1:Nkmx,1:Nmol_max,1:2) integer*8 date1,date2,date_start,date_end integer*8 cputime c integer cput_nbd,cput_init,cput_calc character*(Nchar_mx) datep_file character*(Nchar_mx) command integer code0,codep integer chunk,fpf,fin,nsp character*(Nchar_mx) string,tfile integer Nbin double precision bin(1:Nbin_mx,1:2) integer cdf_indexes(1:Nbin_mx,1:2) double precision cdf(1:Nbin_mx) integer first_bin,last_bin integer current_bin,trigger double precision k_history(0:Nbin_mx) logical sum_lines_contributions logical db double precision epsilon_lr c MPICH c label integer strlen character*(Nchar_mx) label label='program kspectrum' c MPICH call MPI_INIT(ierr) if (ierr.NE.MPI_SUCCESS) THEN write(*,*) 'Error from kspectrum program:' 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,np,ierr) if (np.lt.2) then write(*,*) 'Error from kspectrum program:' write(*,*) 'Number of processes is:',np 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 c MPICH if (pindex.eq.0) then c MPICH write(*,*) 'KSPECTRUM Copyright (C) 2008-2015 Vincent 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(*,*) 'Please consider referencing KSPECTRUM'// & ' in your publications' write(*,*) datep_file='./date' call powercut(over,last_i,last_band,dbr,lastf,lastl,pass) call read_options(sdis,custom_db,htv,htsd_co2,htsd_h2o,dr,nis, & line_profile,usl,sl_choice,slwr,slam,slas, & ucia,ciac,ciawr,ciatr,est,pct,rid, & sens_T,sens_P,sens_x, & uaal,uanb,sda,bslal,ltr,dorg,tfl,rl) call read_data(k_mxep,sd_mxep,nu_min,nu_max, & B_mxep,constant_dnu, & nchunks,T_hitemp,i_inf,i_sup,band_inf,band_sup, & nlw_cz,np_lc,np_bl,tr_dist,ngamma,pr_sd,pr_kc) call init(c2,Rgp,over,rid,band_sigma,band_k, & sens_P,sens_T,sens_x, & band_dsigma_dP,band_dsigma_dT,band_dsigma_dx, & band_dk_dP,band_dk_dT,band_dk_dx) ncalc0=0 ncalc2=0 if (htv.eq.2004) then input_file='./data/molparam_hitran2004.txt' else if (htv.eq.2008) then input_file='./data/molparam_hitran2008.txt' else if (htv.eq.2012) then input_file='./data/molparam_hitran2012.txt' else call error(label) write(*,*) 'Bad input parameter:' write(*,*) 'htv=',htv stop endif call read_molparam(input_file, & Hnmol,Hmol_index,Hmol_niso,Hmol_name, & Hiso_index,Hiso_abundance,Hiso_Qref,Hiso_g,Hiso_mass) call modify_iso_index(Hnmol,Hmol_niso,Hiso_index) call read_composition(nis,Hnmol,Hmol_index,Hmol_niso, & Hiso_index,Hmol_name,Hiso_abundance,Hiso_Qref, & Hiso_g,Hiso_mass, & planet_descriptor,m,alt,P,T, & Nmol_ini,mol_index_ini,mol_niso_ini, & iso_index_ini,mol_name_ini,abundances,iso_Qref_ini, & iso_g_ini,iso_mass_ini,x) call look_for_errors(sdis,custom_db,htv,htsd_co2,htsd_h2o, & dr,nis,line_profile,usl,sl_choice,slwr,slam, & slas,ucia,ciac,ciawr,ciatr, & est,pct,over,rid,sens_T,sens_P,sens_x,uaal,uanb,sda,bslal, & ltr,dorg,tfl,rl,nchunks,k_mxep,sd_mxep,nu_min,nu_max, & B_mxep,constant_dnu, & m,T_hitemp,i_inf,i_sup,band_inf,band_sup, & nlw_cz,np_lc,np_bl,tr_dist,ngamma,pr_sd,pr_kc) call band_limits(nu_min,nu_max,T,m,B_mxep,sdis, & Nb,nu_low,nu_hi) c set i_min,i_max,band_min and band_max call set_limits(over,rid,Nb,m,i_inf,i_sup,band_inf,band_sup, & uaal,uanb,last_i,last_band,i_min,i_max,band_min,band_max, & dbr,lastf,lastl,pass,raz) lastfpp=lastf lastlpp=lastl call print_status(i_min,i_min,i_max,m,band_min,band_min, & band_max,Nb,0,nk,0,0,0,0,0,0.0D+0) do i=i_min,i_max c Debug c write(*,*) 'P(',i,')=',P(i) c write(*,*) 'T(',i,')=',T(i) c Debug call result_file_name(i,kfile,kfile_dP,kfile_dT,kfile_dx) if ((est.eq.1).and.((over.eq.1).or. & ((over.eq.0).and.(.not.rid)))) then command='rm -f '//kfile(1:strlen(kfile)) call exec(command) command='rm -f '//kfile_dP(1:strlen(kfile_dP)) call exec(command) command='rm -f '//kfile_dT(1:strlen(kfile_dT)) call exec(command) command='rm -f '//kfile_dx(1:strlen(kfile_dx)) call exec(command) endif command='rm -f ./results/nu.txt' call exec(command) c set partial pressures of each species for the current atmospheric level call partial_pressures(P,m,x,i,Nmol_ini,mol_index_ini, & mol_niso_ini,iso_index_ini,iso_Qref_ini,iso_g_ini, & iso_mass_ini,mol_name_ini,Ps,Nmol,mol_index, & mol_niso,iso_index,iso_Qref,iso_g,iso_mass,mol_name) c set iso_abundances for the current atmospheric level call set_iso_abundance(abundances,Nmol,mol_niso,i, & iso_abundance) c reading of spectroscopic LBL databases call generate_base(T(i),planet_descriptor,Nmol,mol_index, & iso_index,mol_niso,mol_name,custom_db,htv, & T_hitemp,htsd_co2,htsd_h2o,nis, & Nb,nu_low,nu_hi,over,rid, & nfiles,filename,Tref,fformat,Npass,nl_lp,large_db) call set_Nk_total(raz,i,pass,uanb,Nb,band_inf, & band_min,band_sup,Nk_total,Nk_remain) call set_Sref_frac(sda,P(i),Sref_frac) 555 continue ! multipass loop write(*,48) 'FOR LEVEL',i,' [',i_min,' -',i_max, & ' ], PASS NUMBER:',pass,' /',Npass lastfpp=lastf lastlpp=lastl call calc_Nk_remain(i,pass,Npass,band_min,band_max, & Nk_total,Nk_remain) c Debug c write(*,*) 'calling read_line_parameters' c Debug call read_line_parameters(dr,Nmol,mol_index,mol_niso, & iso_index,nfiles,filename,fformat,pct,rl,pr_kc, & Hiso_abundance,iso_abundance, & Nlines,dbr,lastf,lastl,raz) c Debug c write(*,*) 'read_line_parameters OK' c Debug call write_backup(i,band_min-1,dbr,lastfpp,lastlpp,pass) call get_date(datep_file,date_start) Nk_finished=0 do band=band_min,band_max c Line classification c Debug c write(*,*) 'calling constantk' c Debug call constantk(nu_low(band),nu_hi(band),Nlines,Tref,T,P, & Ps,i,Nmol,mol_index,mol_niso, & iso_index,iso_abundance,iso_mass, & c2,Rgp,iso_g,k_mxep, & line_profile,usl,sl_choice,slwr,slam,slas, & sens_T,sens_P,sens_x, & ltr,dorg,tr_dist,ngamma, & c,sigmac,kc, & dsigmac_dP,dsigmac_dT,dsigmac_dx, & dkc_dP,dkc_dT,dkc_dx, & Nlc,Nlnc,ncl_indexes) c Debug c write(*,*) 'constantk OK' c Dzebug ncalc0=ncalc0+Nlc+Nlnc ncalc2=ncalc2+Nlnc c Narrowband discretization c Debug c write(*,*) 'calling max_error' c Debug call max_error(pass,Nb,band,nu_low, & nu_hi,Sref_frac,Tref,T,P,Ps,m,i,Nmol,mol_index, & mol_niso,iso_index,iso_mass,c2,Rgp,iso_Qref,iso_g, & k_mxep,sd_mxep,pct,sda,constant_dnu, & bslal,nlw_cz,np_lc,np_bl, & sens_P,sens_T,sens_x,pr_sd,large_db, & nk,nu,band_sigma,band_k, & band_dsigma_dP,band_dsigma_dT,band_dsigma_dx, & band_dk_dP,band_dk_dT,band_dk_dx) c Debug c write(*,*) 'max_error OK' c Debug c If truncation over line profiles is required, and if the contribution c of constant lines has to be set to zero, then contribution of constant c lines (sigmac and kc) are NOT added to the result spectrum (goto 412) if ((ltr).and.(tfl)) then goto 412 endif c In all other cases, contribution of constant lines has to be taken into account do ik=1,nk band_sigma(ik)=band_sigma(ik)+sigmac band_k(ik)=band_k(ik)+kc c ------------------------------------------------------------------- c sensitivities if (sens_P) then band_dsigma_dP(ik)=band_dsigma_dP(ik)+dsigmac_dP band_dk_dP(ik)=band_dk_dP(ik)+dkc_dP endif if (sens_T) then band_dsigma_dT(ik)=band_dsigma_dT(ik)+dsigmac_dT band_dk_dT(ik)=band_dk_dT(ik)+dkc_dT endif if (sens_x) then do mol=1,Nmol band_dsigma_dx(ik,mol)=band_dsigma_dx(ik,mol) & +dsigmac_dx(mol) band_dk_dx(ik,mol)=band_dk_dx(ik,mol) & +dkc_dx(mol) enddo ! mol endif c ------------------------------------------------------------------- enddo ! ik 412 continue c Compute physical data needed for computation of k (constant for each narrowband) call precalc(band,sens_P,sens_T,sens_x, & Nlnc,Nmol,i,ncl_indexes,mol_index, & iso_index,mol_niso,Rgp,P,Ps,T,Tref, & iso_abundance,iso_g, & iso_mass,c2, & sorted_index,nuc,dnuc_dP, & gamma_L,dgammaL_dP,dgammaL_dT,dgammaL_dx, & gamma_D,dgammaD_dT, & S,dS_dT, & densities,drho_dP,drho_dT,drho_dx) c Line intensity distribution if (Nlnc.gt.0) then call intensity_distribution(Nlnc,S, & Nbin,bin,cdf_indexes, & cdf,first_bin,last_bin) endif c Debug c write(*,*) 'first_bin=',first_bin c write(*,*) 'last_bin=',last_bin c stop c Debug c Set computation load for each chunk call set_chunks(nchunks,np,nk,ik1,ik2,Nkrun) c write(*,*) 'Now sending data to all child processes' c broadcast data to all child processes c double precision data code0=1 do proc=1,np-1 call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) enddo call MPI_BCAST(i,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(usl,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sl_choice,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(slwr,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(slam,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(slas,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ucia,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(line_profile,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ciac,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ciawr,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ciatr,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sens_T,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sens_P,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sens_x,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(rl,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nmol,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nlnc,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nbin,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(first_bin,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(last_bin,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(cdf_indexes,2*Nbin_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(index,Nline_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sorted_index,Nline_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ncl_indexes,Nline_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(mol_index,Nmol_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(k_mxep,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Ps,Nmol_max,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nuc,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(S,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(densities,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(T,Nmax,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(gamma_L,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(gamma_D,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(bin,2*Nbin_mx,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) if ((sens_P).or.(sens_T).or.(sens_x)) then call MPI_BCAST(dnuc_dP,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dgammaL_dP,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dgammaL_dT,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dgammaL_dx,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dgammaD_dT,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dS_dT,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(drho_dP,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(drho_dT,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(drho_dx,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) endif do proc=1,np-1 free(proc)=1 enddo fin=0 chunk=0 write(*,42) 'Now computing k in narrowband: ',band, & ' [',band_min,' -',band_max,' ]' do while (fin.eq.0) 444 continue call find_free_process(free,np-1,fpf,proc) if (fpf.eq.1) then ! a free child process found chunk=chunk+1 if ((chunk.le.nchunks).and.(Nkrun(chunk).gt.0)) & then do ik=1,Nkrun(chunk) nu_tmp(ik)=nu(ik1(chunk)-1+ik) enddo ch_idx(proc)=chunk code0=2 call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) call MPI_SEND(Nkrun(chunk),1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) call MPI_SEND(nu_tmp,Nkrun(chunk), & MPI_DOUBLE_PRECISION,proc,1, & MPI_COMM_WORLD,ierr) free(proc)=0 ! child process index 'proc' is busy else ! chunk>nchunks or Nkrun(chunk)=0 free(proc)=2 ! means child process index 'proc' has been stopped endif goto 444 else ! no free child process found call stopped_processes(free,np-1,nsp) if (nsp.eq.np-1) then fin=1 goto 111 endif c wait for results call MPI_RECV(proc,1,MPI_INTEGER & ,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD, & stat,ierr) call MPI_RECV(res,Nkmx*2,MPI_DOUBLE_PRECISION & ,proc,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) c ------------------------------------------------------------------- c sensitivities if (sens_P) then call MPI_RECV(dres_dP,Nkmx*2, & MPI_DOUBLE_PRECISION, & proc,MPI_ANY_TAG, & MPI_COMM_WORLD,stat,ierr) endif if (sens_T) then call MPI_RECV(dres_dT,Nkmx*2, & MPI_DOUBLE_PRECISION, & proc,MPI_ANY_TAG, & MPI_COMM_WORLD,stat,ierr) endif if (sens_x) then call MPI_RECV(dres_dx,Nkmx*Nmol_max*2, & MPI_DOUBLE_PRECISION, & proc,MPI_ANY_TAG, & MPI_COMM_WORLD,stat,ierr) endif c ------------------------------------------------------------------- do ik=ik1(ch_idx(proc)),ik2(ch_idx(proc)) c Debug if (res(ik-ik1(ch_idx(proc))+1,2).lt.0.0D+0) & then write(*,*) 'Error from ', & label(1:strlen(label)),' :' write(*,*) 'res(',ik-ik1(ch_idx(proc))+1, & ',2)=', & res(ik-ik1(ch_idx(proc))+1,2) stop endif c Debug band_sigma(ik)=band_sigma(ik) & +res(ik-ik1(ch_idx(proc))+1,1) band_k(ik)=band_k(ik) & +res(ik-ik1(ch_idx(proc))+1,2) c ------------------------------------------------------------------- c sensitivities if (sens_P) then band_dsigma_dP(ik)=band_dsigma_dP(ik) & +dres_dP(ik-ik1(ch_idx(proc))+1,1) band_dk_dP(ik)=band_dk_dP(ik) & +dres_dP(ik-ik1(ch_idx(proc))+1,2) endif if (sens_T) then band_dsigma_dT(ik)=band_dsigma_dT(ik) & +dres_dT(ik-ik1(ch_idx(proc))+1,1) band_dk_dT(ik)=band_dk_dT(ik) & +dres_dT(ik-ik1(ch_idx(proc))+1,2) endif if (sens_x) then do mol=1,Nmol band_dsigma_dx(ik,mol)= & band_dsigma_dx(ik,mol) & +dres_dx(ik-ik1(ch_idx(proc))+1, & mol,1) band_dk_dx(ik,mol)= & band_dk_dx(ik,mol) & +dres_dx(ik-ik1(ch_idx(proc))+1, & mol,2) enddo ! mol endif c ------------------------------------------------------------------- enddo ! ik free(proc)=1 ! child process index 'proc' is free again endif ! fpf=0 or 1 111 continue enddo ! while fin=0 c Debug c write(*,*) 'band_k(',nk,')=',band_k(nk) c Debug c write results to disk, only if it is the last pass ! if (dbr.eq.1) then c Debug c write(*,*) 'In:',label(1:strlen(label)) c write(*,*) 'band_dk_dT=',band_dk_dT(1) c write(*,*) 'band_dk_dP=',band_dk_dP(1) c write(*,*) 'band_dk_dx=',band_dk_dx(1,1) c Debug call write_bandk(Nb,band, & nu_low,nu_hi,m,i,nk,nu, & Nmol,P,T,x,mol_name, & band_sigma,band_k, & sens_P,sens_T,sens_x, & band_dsigma_dP,band_dsigma_dT,band_dsigma_dx, & band_dk_dP,band_dk_dT,band_dk_dx) endif c write temporary results file for each pass call write_bandk_tmp(Nb,band,nu_low,nu_hi,m,i,nk,nu, & Nmol,band_sigma,band_k, & band_dsigma_dP,band_dsigma_dT,band_dsigma_dx, & band_dk_dP,band_dk_dT,band_dk_dx) call write_backup(i,band,dbr,lastfpp,lastlpp,pass) if (pass.eq.1) then Nk_total=Nk_total+nk else Nk_remain=Nk_remain-nk endif datep_file='./date' call get_date(datep_file,date_end) Nk_finished=Nk_finished+nk mtpk=dble(date_end-date_start)/dble(Nk_finished) call print_status(i,i_min,i_max,m,band,band_min,band_max, & Nb,nk,nk,pass,Npass,nl_lp,Nk_total,Nk_remain,mtpk) enddo ! band if (dbr.eq.0) then ! for a new pass pass=pass+1 raz=0 if (uanb.eq.1) then band_min=1 else band_min=band_inf endif call write_backup(i,band_min-1,dbr,lastf,lastl,pass) goto 555 endif c For a new level: pass=1 dbr=0 lastf=1 lastl=0 raz=0 if (uanb.eq.1) then band_min=1 else band_min=band_inf endif c Debug c write(*,*) 'setting band_min=',band_min c Debug c cleaning temporary files that are no longer needed tfile='./optimizations/k*_tmp' call rmfile(tfile) enddo ! i c stopping all child processes do proc=1,np-1 code0=0 call MPI_SEND(code0,1,MPI_INTEGER & ,proc,1,MPI_COMM_WORLD,ierr) call MPI_SEND(pct,1,MPI_LOGICAL & ,proc,1,MPI_COMM_WORLD,ierr) enddo call print_status(i,i_min,i_max,m,band,band_min,band_max, & Nb,nk,nk,Npass,Npass,nl_lp, & Nk_total,Nk_remain,mtpk) call write_info(k_mxep,sd_mxep,ncalc0,ncalc2) call write_over(1) write(*,*) '...done' call clean() else ! every other process string='date' call stringp_file_name(pindex,string,datep_file) cputime=0 333 continue c wait for "codep" from root process call MPI_RECV(codep,1,MPI_INTEGER & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) c codep=0 means exit if (codep.eq.0) then call MPI_RECV(pct,1,MPI_LOGICAL & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) goto 666 c codep=1 means data reception for k computation else if (codep.eq.1) then c receive integer data call MPI_BCAST(i,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(usl,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sl_choice,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(slwr,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(slam,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(slas,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ucia,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(line_profile,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ciac,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ciawr,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ciatr,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sens_T,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sens_P,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(sens_x,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(rl,1,MPI_LOGICAL & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Nmol,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(Nlnc,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(Nbin,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(first_bin,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(last_bin,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(cdf_indexes,2*Nbin_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(index,Nline_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(sorted_index,Nline_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(ncl_indexes,Nline_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(mol_index,Nmol_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) c receive double precision data call MPI_BCAST(k_mxep,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(Ps,Nmol_max,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(nuc,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(S,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(densities,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(T,Nmax,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(gamma_L,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(gamma_D,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(bin,2*Nbin_mx,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) if ((sens_P).or.(sens_T).or.(sens_x)) then call MPI_BCAST(dnuc_dP,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(dgammaL_dP,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(dgammaL_dT,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(dgammaL_dx,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(dgammaD_dT,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(dS_dT,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(drho_dP,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(drho_dT,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(drho_dx,Nlnc,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) endif c codep=2 means a computation of k over a spectral interval must be run else if (codep.eq.2) then call get_date(datep_file,date1) c wait for number of computation and corresponding "nu_tmp" array call MPI_RECV(Nk_p,1,MPI_INTEGER & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) call MPI_RECV(nu_tmp,Nk_p,MPI_DOUBLE_PRECISION & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) do ik=1,Nk_p res(ik,1)=0.0D+0 res(ik,2)=0.0D+0 dres_dP(ik,1)=0.0D+0 dres_dP(ik,2)=0.0D+0 dres_dT(ik,1)=0.0D+0 dres_dT(ik,2)=0.0D+0 do mol=1,Nmol dres_dx(ik,mol,1)=0.0D+0 dres_dx(ik,mol,2)=0.0D+0 enddo if (Nlnc.eq.0) then goto 420 endif k_history(0)=0.0D+0 current_bin=1 trigger=cdf_indexes(current_bin,2) sum_lines_contributions=.true. l=0 if (rl) then epsilon_lr=pr_kc ! value of epsilon when line rejection (for computation of k) is active else epsilon_lr=0.0D+0 ! 0.0D+0: never stop the computation until last line is reached endif c Debug c if (ik.eq.1) then c write(*,*) 'ik=',ik c write(*,*) 'trigger=',trigger c endif c Debug do while (sum_lines_contributions) l=l+1 c compute the contribution of each transition, only if line profile truncation is not c required, or if line profile truncation is required, and distance to line center is c lower than line truncation limit distance. call indexes(Nmol,mol_index,iso_index,mol_niso, & sorted_index(l),1,molec,isotop,isof) call calc_k_opt(nu_tmp(ik),nuc(l),dnuc_dP(l), & P(i),Ps(molec),T(i),sorted_index(l), & gamma_L(l),dgammaL_dP(l), & dgammaL_dT(l),dgammaL_dx(l), & gamma_D(l),dgammaD_dT(l), & S(l),dS_dT(l), & densities(l),drho_dP(l),drho_dT(l),drho_dx(l), & line_profile,usl,sl_choice,slwr,slam,slas, & sens_T,sens_P,sens_x, & sigma,k,dsigma_dP,dsigma_dT,dsigma_dx, & dk_dP,dk_dT,dk_dx) res(ik,1)=res(ik,1)+sigma ! only allowed transitions res(ik,2)=res(ik,2)+k ! only allowed transitions c ------------------------------------------------------------------- c sensitivities (allowed transitions only) if (sens_P) then dres_dP(ik,1)=dres_dP(ik,1)+dsigma_dP dres_dP(ik,2)=dres_dP(ik,2)+dk_dP endif if (sens_T) then dres_dT(ik,1)=dres_dT(ik,1)+dsigma_dT dres_dT(ik,2)=dres_dT(ik,2)+dk_dT endif if (sens_x) then call indexes(Nmol,mol_index,iso_index,mol_niso, & sorted_index(l),1,molec,isotop,isof) dres_dx(ik,molec,1)=dres_dx(ik,molec,1) & +dsigma_dx dres_dx(ik,molec,2)=dres_dx(ik,molec,2) & +dk_dx endif c ------------------------------------------------------------------- if (l.eq.trigger) then call line_intervals(db,l,Nbin,bin,cdf_indexes, & res(ik,2),k_mxep,epsilon_lr, & current_bin,trigger, & k_history,sum_lines_contributions) c Debug c if (ik.eq.1) then c write(*,*) 'l=',l,' sum_lines_contributions=', c & sum_lines_contributions, c & ' trigger=',trigger c endif c Debug endif enddo ! while (sum_lines_contributions) 420 continue c Debug c if (ik.eq.1) then c write(*,*) 'sum_lines_contributions=', c & sum_lines_contributions c write(*,*) 'l=',l c endif c Debug c Collision-induced absorption (for CO2 only) co2f=0 do mol=1,Nmol if (mol_index(mol).eq.2) then co2f=1 goto 421 endif enddo 421 continue if (co2f.eq.1) then call CIA(T(i),nu_tmp(ik),Nmol,P(i),Ps,mol_index, & ucia,ciac,ciawr,ciatr, & sens_P,sens_T,sens_x, & kCIA,dkCIA_dP,dkCIA_dT,dkCIA_dxco2) res(ik,2)=res(ik,2)+kCIA ! [m¯¹] allowed transitions + CIA transitions c ------------------------------------------------------------------- c sensitivities (CIA transitions) if (sens_P) then dres_dP(ik,2)=dres_dP(ik,2)+dkCIA_dP endif if (sens_T) then dres_dT(ik,2)=dres_dT(ik,2)+dkCIA_dT endif if (sens_x) then do mol=1,Nmol if (mol_index(mol).eq.2) then dres_dx(ik,mol,2)=dres_dx(ik,mol,2) & +dkCIA_dxCO2 endif enddo ! mol endif c ------------------------------------------------------------------- endif enddo ! ik c send results to root process call MPI_SEND(pindex,1,MPI_INTEGER & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(res,Nkmx*2,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) c ------------------------------------------------------------------- c sensitivities if (sens_P) then call MPI_SEND(dres_dP,Nkmx*2,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) endif if (sens_T) then call MPI_SEND(dres_dT,Nkmx*2,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) endif if (sens_x) then call MPI_SEND(dres_dx,Nkmx*Nmol_max*2, & MPI_DOUBLE_PRECISION, & 0,1,MPI_COMM_WORLD,ierr) endif c ------------------------------------------------------------------- call get_date(datep_file,date2) cputime=cputime+(date2-date1) c codep=3 means data reception for spectral discretization else if (codep.eq.3) then c receive integer data call MPI_BCAST(sda,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(bslal,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(Nmol,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(i,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(mol_index,Nmol_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(mol_niso,Nmol_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(iso_index,Nmol_max*Niso_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(iso_g,Nmol_max*Niso_max,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(nacc,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(nxtab,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(nk_func,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(nx,Nytab_mx,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(kindex,Nytab_mx*2,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(nlw_cz,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(np_lc,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(np_bl,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,stat,ierr) c receive double precision data call MPI_BCAST(constant_dnu,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(c2,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(Sref_frac,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(iso_mass,Nmol_max*Niso_max, & MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(Tref,Nmol_max,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(T,Nmax,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(P,Nmax,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(Ps,Nmol_max,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(xtab,Nxtab_mx*3,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(xacc,Nacc_mx*2,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(ktab,Nktab_mx*3,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(yvalues,Nytab_mx,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) call MPI_BCAST(pr_sd,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,stat,ierr) c codep=4 means computation for spectral discretization else if (codep.eq.4) then call get_date(datep_file,date1) c receive data call MPI_RECV(numin,1,MPI_DOUBLE_PRECISION & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) call MPI_RECV(numax,1,MPI_DOUBLE_PRECISION & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) c compute call max_error_1cm_alrl(pindex,numin,numax,c2, & sda,constant_dnu,bslal, & nlw_cz,np_lc,np_bl,Sref_frac,Nmol,iso_g, & mol_index,iso_index,mol_niso,Ps,Tref,T,P,i,iso_mass, & nacc,nxtab,nk_func,nx,kindex, & xtab,xacc,ktab,yvalues,pr_sd, & nk_tmp,nu_tmp) c send results call MPI_SEND(pindex,1,MPI_INTEGER & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(nk_tmp,1,MPI_INTEGER & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(nu_tmp,Nkmx,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) call get_date(datep_file,date2) cputime=cputime+(date2-date1) c codep=5 means data reception for voigt tabulation else if (codep.eq.5) then c receive integer data call MPI_BCAST(Nit_max,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(ndiv,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) c receive double precision data call MPI_BCAST(eps_max,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dx0_ori,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(pp,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(dp_max,1,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) c codep=6 means computation for voigt tabulation else if (codep.eq.6) then call get_date(datep_file,date1) c receive data call MPI_RECV(y,1,MPI_DOUBLE_PRECISION & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) c compute dx0=dx0_ori call voigt_tabulation_1y(y,eps_max,dx0,Nit_max,pp,ndiv, & dp_max,n1,n2,xres) c send results call MPI_SEND(pindex,1,MPI_INTEGER & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(n1,1,MPI_INTEGER & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(n2,1,MPI_INTEGER & ,0,1,MPI_COMM_WORLD,ierr) call MPI_SEND(xres,Nppc_max*3,MPI_DOUBLE_PRECISION & ,0,1,MPI_COMM_WORLD,ierr) call get_date(datep_file,date2) cputime=cputime+(date2-date1) c codep=7 means data reception for file parsing else if (codep.eq.7) then c receive integer data call MPI_BCAST(nsint,1,MPI_INTEGER & ,0,MPI_COMM_WORLD,ierr) c receive double precision data call MPI_BCAST(slimits1,Nti_mx*2,MPI_DOUBLE_PRECISION & ,0,MPI_COMM_WORLD,ierr) c codep=8 means computation for file parsing else if (codep.eq.8) then call get_date(datep_file,date1) call MPI_RECV(ifile,1,MPI_INTEGER & ,0,MPI_ANY_TAG,MPI_COMM_WORLD,stat,ierr) call parse_file(ifile,nsint,slimits1) c send results call MPI_SEND(pindex,1,MPI_INTEGER & ,0,1,MPI_COMM_WORLD,ierr) call get_date(datep_file,date2) cputime=cputime+(date2-date1) endif ! codep c wait for a new execution code goto 333 666 continue if (pct) then write(*,43) 'Total computation time for process', pindex, & ' is :' call print_cputime(cputime) endif endif ! pindex c MPICH call MPI_FINALIZE(ierr) c MPICH end