diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e488c13 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*_mpitmp.f +*.o +*.a +xcql3d* \ No newline at end of file diff --git a/ADC/Makefile b/ADC/Makefile new file mode 100755 index 0000000..375bb6a --- /dev/null +++ b/ADC/Makefile @@ -0,0 +1,14 @@ +# Makefile that builts a library lib$(LIB).a from all +# of the Fortran files found in the current directory. +# Usage: make LIB= +# Pearu +LIB=ADC +OBJ=$(patsubst %.f,%.o,$(shell ls ADC/*.f)) +all: lib$(LIB).a +$(OBJ): + $(FC) -c $(FFLAGS) $(patsubst ADC/%.o,ADC/%.f,ADC/$(@F)) -o $@ +lib$(LIB).a: $(OBJ) + $(AR) rus ADC/lib$(LIB).a $? +clean: + rm -f ADC/*.o + rm -f ADC/*.a \ No newline at end of file diff --git a/ainpla.f b/ainpla.f index 79f9aa1..d45f5fa 100644 --- a/ainpla.f +++ b/ainpla.f @@ -1836,101 +1836,102 @@ subroutine impurity_update CMPIINSERT_INCLUDE include 'param.h' include 'comm.h' + external adcdo save dMpellet_sum - real*8 ZDISTR(100) ! local, for subr. ADCDO() + real*8 ZDISTR(100) ! local, for subr. ADCDO() - if((imp_depos_method.ne.'disabled') .and. (kelecg.eq.1))then + if((imp_depos_method.ne.'disabled') .and. (kelecg.eq.1))then !YuP[2020-06-24] Changed (gamafac.eq."hesslow") to (imp_depos_method.ne.'disabled') ! [a more general logic] - !YuP[2019-09-16] - ! After temp() and reden() profiles are updated, - ! find profile of impurity produced by pellet: - if(imp_depos_method.eq.'pellet' - & .and. timet.ge.pellet_tstart) then - if(dMpellet_sum.lt.pellet_M0)then ! Not all pellet yet gone - kopt=1 !propagate pellet; calc. dMpellet_dvol_sum(1:lrz) - Rpellet= pellet_Rstart -pellet_V*(timet-pellet_tstart) - if(Rpellet.lt.rmag)then - ! ablation rate is determined by ne and Te - ! which are already affected by pellet deposition - temp_wk(1:lrz)=temp(kelecg,1:lrz) - reden_wk(1:lrz)=reden(kelecg,1:lrz) - else ! Rpellet>rmag - !use values before pellet affected the temp. and dens. - temp_wk(1:lrz)=temp_t0(kelecg,1:lrz) - reden_wk(1:lrz)=reden_t0(kelecg,1:lrz) - endif - call set_get_pellet(kopt, timet, lrz, rmag, - & rpcon(1:lrz), rmcon(1:lrz), dvol(1:lrz), - & temp_wk(1:lrz), reden_wk(1:lrz), - & pellet_V, pellet_M0, pellet_Rstart, pellet_tstart, - & pellet_rcloud1, pellet_rcloud2, pellet_rp0, - & pellet_pn, pellet_pt, pellet_pm, - & ipellet_method, pellet_fract_rem, - & ipellet_iter_max, pellet_iter_accur, pellet_Cablation, - & Gablation,dMpellet_sum, dMpellet_dvol_sum(1:lrz) ) -!CMPIINSERT_IF_RANK_EQ_0 -! do ll=1,lrz -! WRITE(*,'(a,i4,2e14.7)') -! + 'tdchief/set_get_pellet: n,dMpellet_sum,dMpellet_dvol_sum(ll)', -! + n,dMpellet_sum,dMpellet_dvol_sum(ll) -! enddo -!CMPIINSERT_ENDIF_RANK - - ! Find and save local rho at pellet position (for plots) - rloc= Rpellet ! Local R at pellet position - zloc= zmag ! Assumed - pellet travels along Z=Zaxis - ppsi=terp2(rloc,zloc,nnr,er,nnz,ez,epsi,epsirr,epsizz, - & epsirz,nnra,0,0) - ! From ppsi, find local rho, based on eqrho() array . - ! Definition of eqrho is set by specification of radcoord. - call terp1(nconteqn,-eqpsi,eqrho,d2eqrho,-ppsi,1,tab,itab) - rho_loc=tab(1)/rhomax ! rhomax is in comm.h - pellet_rho= rho_loc - pellet_Mrem= pellet_M0-dMpellet_sum ! Remaining mass[gram] - !Next: calculate particle density, use convert=Avogadro/atw - ! where atw is the atomic weight (40g/mol for Argon), - ! or simply (Avogadro*proton)/(atw*proton) = 1.0/fmass_imp - ! where fmass_imp=atw*proton (in comm.h) - ! THIS IS WHAT WE NEED FOR cfpcoefn.f (stored in comm.h): - dens_imp_allstates(1:lrz)=dMpellet_dvol_sum(1:lrz)/fmass_imp - ! Note: At given instant t=timet, pellet is at - ! Rpellet(t)= pellet_Rstart -pellet_V*(timet-pellet_tstart) - ! It means that a given flux surface at R=rpcon(lr) - ! is reached by pellet at instant - ! t= pellet_tstart +(pellet_Rstart-rpcon(lr))/pellet_V - ! We use this t as the onset time for the - ! exp(-(t-tstart)/tau) temperature drop, - ! i.e. as the value for - ! temp_expt_tstart(lr) in case of iprote='prb-expt' option. - ! Similarly in case of iproti='prb-expt', - ! and we use same values of temp_expt_tstart, etc. - endif !(dMpellet_sum.lt.pellet_M0) Not all pellet yet gone - ! But what to do when ALL pellet is ablated? - ! For now, we assume, all that pellet material stays within - ! flux surfaces (no radial transport). Only ionization - ! states may change, because of changing Te(time). - endif !(imp_depos_method.eq.'pellet') - - if(imp_depos_method.eq.'instant' - & .and. timet.ge.tstart_imp) then !YuP[2019-12-05] - !Instant deposition of impurity at all surfaces. - ! dens0_imp(0:lrz) profile must be set in namelist - dens_imp_allstates(1:lrz)=dens0_imp(1:lrz) !YuP[2019-12-05] - endif !(imp_depos_method.eq.'instant') !YuP[2019-12-05] - - !--------- NOW Find distribution over charge states. - ! Based on Corona model, and using ADPAK data (from *.ntau files). - ! general INPUT arguments for sub.set_get_ADPAK() : - kopt=1 ! option flag: '1' is to get and - tau_r=adpak_tau_r ![sec] Characteristic time of radial decay of T_e - do ll=1,lrz !ilend - call tdnflxs(ll) !-> get l_,lr_,... - dens_ne=reden(kelecg,lr_) ![cm^-3] - temp_Te= temp(kelecg,lr_) ! Te [keV] at this flux surface - !Get neutral density of D0, needed as an input for ADPAK tables, - ! where it enters through the ratio nD0/ne. - rho=rya(lr_) + !YuP[2019-09-16] + ! After temp() and reden() profiles are updated, + ! find profile of impurity produced by pellet: + if(imp_depos_method.eq.'pellet' + & .and. timet.ge.pellet_tstart) then + if(dMpellet_sum.lt.pellet_M0)then ! Not all pellet yet gone + kopt=1 !propagate pellet; calc. dMpellet_dvol_sum(1:lrz) + Rpellet= pellet_Rstart -pellet_V*(timet-pellet_tstart) + if(Rpellet.lt.rmag)then + ! ablation rate is determined by ne and Te + ! which are already affected by pellet deposition + temp_wk(1:lrz)=temp(kelecg,1:lrz) + reden_wk(1:lrz)=reden(kelecg,1:lrz) + else ! Rpellet>rmag + !use values before pellet affected the temp. and dens. + temp_wk(1:lrz)=temp_t0(kelecg,1:lrz) + reden_wk(1:lrz)=reden_t0(kelecg,1:lrz) + endif + call set_get_pellet(kopt, timet, lrz, rmag, + & rpcon(1:lrz), rmcon(1:lrz), dvol(1:lrz), + & temp_wk(1:lrz), reden_wk(1:lrz), + & pellet_V, pellet_M0, pellet_Rstart, pellet_tstart, + & pellet_rcloud1, pellet_rcloud2, pellet_rp0, + & pellet_pn, pellet_pt, pellet_pm, + & ipellet_method, pellet_fract_rem, + & ipellet_iter_max, pellet_iter_accur, pellet_Cablation, + & Gablation,dMpellet_sum, dMpellet_dvol_sum(1:lrz) ) +!CMPIINSERT_IF_RANK_EQ_0 +! do ll=1,lrz +! WRITE(*,'(a,i4,2e14.7)') +! + 'tdchief/set_get_pellet: n,dMpellet_sum,dMpellet_dvol_sum(ll)', +! + n,dMpellet_sum,dMpellet_dvol_sum(ll) +! enddo +!CMPIINSERT_ENDIF_RANK + + ! Find and save local rho at pellet position (for plots) + rloc= Rpellet ! Local R at pellet position + zloc= zmag ! Assumed - pellet travels along Z=Zaxis + ppsi=terp2(rloc,zloc,nnr,er,nnz,ez,epsi,epsirr,epsizz, + & epsirz,nnra,0,0) + ! From ppsi, find local rho, based on eqrho() array . + ! Definition of eqrho is set by specification of radcoord. + call terp1(nconteqn,-eqpsi,eqrho,d2eqrho,-ppsi,1,tab,itab) + rho_loc=tab(1)/rhomax ! rhomax is in comm.h + pellet_rho= rho_loc + pellet_Mrem= pellet_M0-dMpellet_sum ! Remaining mass[gram] + !Next: calculate particle density, use convert=Avogadro/atw + ! where atw is the atomic weight (40g/mol for Argon), + ! or simply (Avogadro*proton)/(atw*proton) = 1.0/fmass_imp + ! where fmass_imp=atw*proton (in comm.h) + ! THIS IS WHAT WE NEED FOR cfpcoefn.f (stored in comm.h): + dens_imp_allstates(1:lrz)=dMpellet_dvol_sum(1:lrz)/fmass_imp + ! Note: At given instant t=timet, pellet is at + ! Rpellet(t)= pellet_Rstart -pellet_V*(timet-pellet_tstart) + ! It means that a given flux surface at R=rpcon(lr) + ! is reached by pellet at instant + ! t= pellet_tstart +(pellet_Rstart-rpcon(lr))/pellet_V + ! We use this t as the onset time for the + ! exp(-(t-tstart)/tau) temperature drop, + ! i.e. as the value for + ! temp_expt_tstart(lr) in case of iprote='prb-expt' option. + ! Similarly in case of iproti='prb-expt', + ! and we use same values of temp_expt_tstart, etc. + endif !(dMpellet_sum.lt.pellet_M0) Not all pellet yet gone + ! But what to do when ALL pellet is ablated? + ! For now, we assume, all that pellet material stays within + ! flux surfaces (no radial transport). Only ionization + ! states may change, because of changing Te(time). + endif !(imp_depos_method.eq.'pellet') + + if(imp_depos_method.eq.'instant' + & .and. timet.ge.tstart_imp) then !YuP[2019-12-05] + !Instant deposition of impurity at all surfaces. + ! dens0_imp(0:lrz) profile must be set in namelist + dens_imp_allstates(1:lrz)=dens0_imp(1:lrz) !YuP[2019-12-05] + endif !(imp_depos_method.eq.'instant') !YuP[2019-12-05] + + !--------- NOW Find distribution over charge states. + ! Based on Corona model, and using ADPAK data (from *.ntau files). + ! general INPUT arguments for sub.set_get_ADPAK() : + kopt=1 ! option flag: '1' is to get and + tau_r=adpak_tau_r ![sec] Characteristic time of radial decay of T_e + do ll=1,lrz !ilend + call tdnflxs(ll) !-> get l_,lr_,... + dens_ne=reden(kelecg,lr_) ![cm^-3] + temp_Te= temp(kelecg,lr_) ! Te [keV] at this flux surface + !Get neutral density of D0, needed as an input for ADPAK tables, + ! where it enters through the ratio nD0/ne. + rho=rya(lr_) !YuP[2020-11-15] Moved call to get_dens_nD0_ADPAK() outside of if(). ! It is not actually specific to ADPAK. ! It simply sets dens_nD0= dens_nD0_b*exp((rho-1.d0)*radmin/dens_nD0_l) @@ -1939,56 +1940,56 @@ subroutine impurity_update ![2020-11-15] This change affects results, but not strongly. ! Before this change the value of dens_nD0 was not defined ! (initialized to a neg. value on PC, or small pos. value on Cray) - if(adpak.eq.'enabled')then - call set_get_ADPAK(kopt,imp_type, - & temp_Te,dens_nD0,dens_ne,tau_r, - & z1av,z2av) !-> OUT: and for given imp_type - !Now, from knowledge of and , find distribution fz(Z) - ! over charge states: - iZatom=INT(bnumb_imp(nstates)) !should be integer to this subr. - !Could simply use iZatom=nstates - call get_distr_charge_states(iZatom,z1av,z2av, fz) - ! INPUT: iZatom,z1av,z2av OUTPUT: fz(0:iZatom) - !Note: fz has normalization SUM(fz(0:iZatom))=1, - ! and SUM(fz*Z)= (with some accuracy). - else ! adpak.ne.'enabled' - ! Alternatively (use it for W; which is imp_type=8): - INUCZ=INT(bnumb_imp(nstates)) !atom nuclear charge - ZTE= temp_Te*1d3 ! ZTE is Te in eV - ZNE= dens_ne ! electron density in cm-3 - ZNA= dens_nD0 ! density of hydrogen atoms in cm-3 - kion1=kionm(1) !assume all ions have same T, so pick the 1st - ZTA=temp(kion1,lr_)*1d3/2. - !ZTA= Timp/Mimp +Ta/Ma ~ Ta/Ma. For D: Ti/2 (Ma=2) - !ZTA is relative temperature for charge exchange of ions with - !hydrogen atoms which is Tz/Mz+Th/Mh - CALL ADCDO(INUCZ,ZTE,ZNE,ZTA,ZNA, - & ZS, ZS2, EZRAD, EZLOSS, ZDISTR) - !output: - ! ZS is ; ZS2 is - ! EZRAD is the radiation power W*cm3 - ! EZLOSS is the ionization energy losses, W*cm3 - ! ZDISTR is array for distribution function of excited states with - ! sum(ZDISTR)=1 - ! The ADCDO subroutine is in /ADC/ folder. Need to compile and link - ! all 16 files in that folder. - fz(0:nstates)=ZDISTR(1:nstates+1) - endif ! adpak= enabled or not. - !For a given total density of impurity dens_imp_allstates(lr_) - !we now can find the density for each charge state: - do kstate=0,nstates - dens_imp(kstate,lr_)=dens_imp_allstates(lr_)*fz(kstate) !cm^-3 - enddo - !dens_imp_allstates(lr_) is the total - !(over all charge states, incl. Z=0) - !density of impurity at given lr_, at given time step. - ! This is just the ablated material from pellet, - ! before ionization process occured. - !Note: sum(dens_imp(0:nstates,lr_))= dens_imp_allstates(lr_) - enddo ! do ll - call profiles ! update reden(kelecg) based on dens_imp() - !(very small effect from this update of reden) - endif !((imp_depos_method.ne.'disabled') .and. kelecg.eq.1) - + if(adpak.eq.'enabled')then + call set_get_ADPAK(kopt,imp_type, + & temp_Te,dens_nD0,dens_ne,tau_r, + & z1av,z2av) !-> OUT: and for given imp_type + !Now, from knowledge of and , find distribution fz(Z) + ! over charge states: + iZatom=INT(bnumb_imp(nstates)) !should be integer to this subr. + !Could simply use iZatom=nstates + call get_distr_charge_states(iZatom,z1av,z2av, fz) + ! INPUT: iZatom,z1av,z2av OUTPUT: fz(0:iZatom) + !Note: fz has normalization SUM(fz(0:iZatom))=1, + ! and SUM(fz*Z)= (with some accuracy). + else ! adpak.ne.'enabled' + ! Alternatively (use it for W; which is imp_type=8): + INUCZ=INT(bnumb_imp(nstates)) !atom nuclear charge + ZTE= temp_Te*1d3 ! ZTE is Te in eV + ZNE= dens_ne ! electron density in cm-3 + ZNA= dens_nD0 ! density of hydrogen atoms in cm-3 + kion1=kionm(1) !assume all ions have same T, so pick the 1st + ZTA=temp(kion1,lr_)*1d3/2. + !ZTA= Timp/Mimp +Ta/Ma ~ Ta/Ma. For D: Ti/2 (Ma=2) + !ZTA is relative temperature for charge exchange of ions with + !hydrogen atoms which is Tz/Mz+Th/Mh + CALL ADCDO(INUCZ,ZTE,ZNE,ZTA,ZNA, + & ZS, ZS2, EZRAD, EZLOSS, ZDISTR) + !output: + ! ZS is ; ZS2 is + ! EZRAD is the radiation power W*cm3 + ! EZLOSS is the ionization energy losses, W*cm3 + ! ZDISTR is array for distribution function of excited states with + ! sum(ZDISTR)=1 + ! The ADCDO subroutine is in /ADC/ folder. Need to compile and link + ! all 16 files in that folder. + fz(0:nstates)=ZDISTR(1:nstates+1) + endif ! adpak= enabled or not. + !For a given total density of impurity dens_imp_allstates(lr_) + !we now can find the density for each charge state: + do kstate=0,nstates + dens_imp(kstate,lr_)=dens_imp_allstates(lr_)*fz(kstate) !cm^-3 + enddo + !dens_imp_allstates(lr_) is the total + !(over all charge states, incl. Z=0) + !density of impurity at given lr_, at given time step. + ! This is just the ablated material from pellet, + ! before ionization process occured. + !Note: sum(dens_imp(0:nstates,lr_))= dens_imp_allstates(lr_) + enddo ! do ll + call profiles ! update reden(kelecg) based on dens_imp() + !(very small effect from this update of reden) + endif !((imp_depos_method.ne.'disabled') .and. kelecg.eq.1) + return end subroutine impurity_update diff --git a/makefile_gfortran_mpi_rhel7 b/makefile_gfortran_mpi_rhel7 new file mode 100644 index 0000000..d39dd72 --- /dev/null +++ b/makefile_gfortran_mpi_rhel7 @@ -0,0 +1,149 @@ +# makefile for CQL3D using gfortran64 for mpi on PC + +# May/06: makefile_mpi.xxx is a system setup by Nikolai Ershov +# for creation of mpi version of cql3d from the serial +# version. Therefore, only one version of cql3d needs +# be maintained. +# Comments (beginning with CMPI) have been inserted in the source. +# These are processed using python code from the ./mpi subdirectory +# to produce an MPI cql3d version. +# The procedure is described in CompX report +# CQL3D_Parallelization_CompX-2006-1.pdf, and has been evolved +# by Yuri Petrov. +# Please follow this procedure in future cql3d modifications. + +SHELL = /bin/sh +NAME = xcql3d_mpi.gfortran64 +COMPILER= mpif77 +BUILDER= gfortran +INCLUDES = advnce.h comm.h frcomm.h frname.h frname_decl.h name.h \ + param.h trans.h wpadvnc.h name_decl.h +SOURCES = $(wildcard ADC/!(*mpitmp*).f) +SOURCES += a_cqlp.f abchief.f achief1.f achiefn.f aclear.f ainalloc.f \ + aindflt.f aindflt1.f aindfpa.f aingeom.f ainitial.f \ + ainpla.f ainplt.f ainpltpa.f \ + ainsetpa.f ainsetva.f ainspec.f ainvnorm.f ampfar.f \ + aminmx.f bavdens.f bavgmax.f baviorbt.f \ + bcast.f bsl.f bsu.f cfpcoefc.f cfpcoefn.f \ + cfpcoefr.f cfpgamma.f cfpleg.f \ + cfpmodbe.f cfpsymt.f coefefad.f coefefld.f coefegad.f \ + coeffpad.f coefload.f coefmidt.f coefmidv.f \ + coefrfad.f coefstup.f coefsyad.f coefwti.f \ + coefwtj.f diag.f diagcfac.f \ + diagdens.f diagdenz.f diagentr.f diagescl.f diaggnde.f \ + diaggnde2.f diagimpd.f diagscal.f diagwrng.f diagxswt.f \ + diagxswx.f dsk_gr.f dskout.f efield.f eflditer.f eqalloc.f \ + eqcoord.f eqelpse.f eqflxavg.f eqfn.f \ + eqfndpsi.f eqfninv.f eqfpsi.f \ + eqindflt.f eqinitl.f eqjac.f eqonovrp.f \ + eqorbit.f eqrhopsi.f eqrhs.f \ + eqtopeol.f equilib.f eqvolpsi.f eqwrng.f \ + esefld.f exlin.f exsweep.f exsweept.f \ + exsweepx.f finit.f firstdrv.f fle.f flxfn.f \ + freya.f freyasou.f frhexdrv.f \ + frinitl.f frinitz.f frnbdep2.f frnfreya.f \ + frplteq.f frset.f frsmooth.f frsplft.f \ + frstup.f frsubs.f frsuppor.f frwrong.f hpalloc0.f \ + ilut.f impavnc0.f impchk.f impnorm.f it3dalloc.f \ + lookup.f losscone.f lossegy.f lossorbm.f \ + losstor.f micfrplt.f micgetr.f \ + micgmbnd.f micgnbnd.f micxinil.f micxinim.f micxinit.f \ + micxiniz.f netcdfrf.f netcdfrw2.f ntdstore.f \ + ntloop.f pack21.f pltcont.f pltcycl.f pltdf.f \ + pltdnz.f pltelec.f pltendn.f \ + pltends.f pltfluxs.f pltfofvv.f pltfvsv.f pltinit.f pltlosc.f \ + pltmain.f pltpower.f pltprppr.f \ + pltrstv.f pltrun.f pltstrml.f pltvec.f pltvectr.f \ + pltvflux.f profaxis.f profiles.f prppr.f \ + prpprctr.f psif.f psifp.f psifppy.f psifpy.f psiinv.f \ + r8lsode.f r8subs.f rdc_multi.f rdc_bplt.f restcon.f resthks.f \ + restvty.f rf.f sigalloc.f siggy.f sigmax.f sigsetup.f \ + sigv5d.f sigfn.f sigie.f sigmaxwl.f sigv.f \ + soucrit.f sounorm.f soup.f soup0.f souplt.f \ + sourc0.f sourcee.f sourcef.f sourceko.f sourcpwr.f \ + synchrad.f tdbootst.f tdboothi.f tdchief.f tddiag.f \ + tddsig.f tdeqdsk.f tdfinterp.f tdinitl.f tdinlegw.f \ + tdinterp.f tdnflxs.f tdnpa.f tdnpadiag.f tdnpa0.f \ + tdnpacxcs.f tdnpalam.f tdnpabscs.f tdoutput.f \ + tdplteq.f tdpltjop.f tdpltmne.f tdpro.f \ + tdreadf.f tdrmshst.f tdsetnpa.f tdsetsxr.f \ + tdstin.f tdsxr.f tdsxr0.f tdsxray.f tdsxrplt.f \ + tdtloop.f tdtoarad.f \ + tdtoaray.f tdtraloc.f tdtransp.f tdtranspn.f tdtravct.f \ + tdtrchk.f tdtrchkd.f tdtrcon.f \ + tdtrdfus.f tdtrfcop.f tdtrflg.f tdtrflx.f \ + tdtrmuy.f tdtrrsou.f tdtrrtov.f \ + tdtrrtov2.f tdtrsavf.f tdtrsym.f tdtrvint.f \ + tdtrvsou.f tdtrvtor.f tdtrvtor2.f \ + tdtrvtor3.f tdtrwtl.f tdtry.f tdtscinp.f tdtscout.f \ + tdwrng.f tdwritef.f tdxin13d.f \ + tdxin23d.f tdxin33d.f tdxinitl.f urfalloc.f \ + urfavg.f urfb0.f urfbes.f \ + urfbplt.f urfchief.f urfdamp0.f urfdamp1.f \ + urfdamp2.f urfdampa.f urfdout.f urfedge.f \ + urffflx.f urfindfl.f urfinitl.f urfmidv.f \ + urfpack.f urfpackm.f urfrays.f urfread.f \ + urfread_.f urfsetup.f urfwrite.f urfwrite_.f \ + urfwrong.f urfwr0.f urfwr0c.f \ + vlf.f vlfalloc.f vlfbplt.f vlfsetup.f vlh.f vlhbplt.f vlhd.f \ + wpalloc.f wparsou.f wpavg.f wpbdry.f \ + wpcheck.f wpchgdy.f wpcthta.f wpelecf.f \ + wpinitl.f wploweq.f wpsavf.f wptrafx.f wptramu.f \ + wptrmuy.f wpvptb.f wpwrng.f wpmshchk.f \ + zblock.f zcunix.f zfreya.f znonsym.f +PYTHON = python3 +OBJECTS = $(SOURCES:.f=.o) +LOCATION = -L$(NETCDF_DIR)/lib -L$(PGPLOT_DIR)/lib -L$(PWD)/ADC +LIBS= -lX11 -lnetcdff -lnetcdf -lpgplot -lADC +INCLUDE=/usr/include +#DEBUG = -g -Mbounds +DEBUG = -g +OPTIMIZE = -O +LISTING = -Mlist +CSPECIAL = -Wno-padded +FFLAGS = $(CSPECIAL) $(OPTIMIZE) +#SPECIAL = -byteswapio -Mbackslash +LDSPECIAL = -Wl,-noinhibit-exec,-rpath,$(PGPLOT_DIR)/lib +#COMPILE = mpif77 -c $(CSPECIAL) -c $(DEBUG) -I $(INCLUDE) # or use $(OPTIMIZE) +COMPILE = $(COMPILER) -c $(FFLAGS) -I $(INCLUDE) # or use $(DEBUG) +LOAD = $(COMPILER) -o $(NAME) $(LDSPECIAL) # Remove -m for optimize +PROTECT = chmod 755 +DELETE = rm -f + +# The following gives suffixes to be used in checking for suffix rules. +# Written without dependencies, it may be useful to turn of such checking? +.SUFFIXES: + +$(NAME): $(OBJECTS) + $(LOAD) $(OBJECTS) $(LOCATION) $(LIBS) + $(PROTECT) $(NAME) + +# Following use of pattern matching works; +# it is based on statements in Sect. 10.7 of gmake Manual +# by Stallman and McGrath. +# However, the two rules following this seem clearer to me (BobH). +#%.o: %.f $(INCLUDES) +# $(COMPILE) $< -o $@ + +include deps #Introduced by John Wright (JCW) +$(SOURCES): $(INCLUDES) + +# see Sect. 4.10.1, Static Pattern Rules. +#$(OBJECTS): %.o: %.f $(INCLUDES) +%.o: %.f # $(INCLUDES) +#Changes introduced by JCW so obtain saved copy of mpi modified sources +# mpi/doparallel.py $< mpitmp.f mpi/mpins_par.f +# $(COMPILE) mpitmp.f -o $@ + export COMPILER FFLAGS + make -f ADC/Makefile + $(PYTHON) mpi/doparallel.py $< $*_mpitmp.f mpi/mpins_par.f + $(COMPILE) $*_mpitmp.f -o $@ + +rebuild: + $(COMPILE) $(SOURCES) + $(LOAD) $(OBJECTS) $(LOCATION) $(LIBS) + +clean: + make -f ADC/Makefile clean + $(DELETE) $(OBJECTS) *.lst *_mpitmp.f + diff --git a/makefile_mpi.gfortran64_rhel b/makefile_mpi.gfortran64_rhel new file mode 100644 index 0000000..d39dd72 --- /dev/null +++ b/makefile_mpi.gfortran64_rhel @@ -0,0 +1,149 @@ +# makefile for CQL3D using gfortran64 for mpi on PC + +# May/06: makefile_mpi.xxx is a system setup by Nikolai Ershov +# for creation of mpi version of cql3d from the serial +# version. Therefore, only one version of cql3d needs +# be maintained. +# Comments (beginning with CMPI) have been inserted in the source. +# These are processed using python code from the ./mpi subdirectory +# to produce an MPI cql3d version. +# The procedure is described in CompX report +# CQL3D_Parallelization_CompX-2006-1.pdf, and has been evolved +# by Yuri Petrov. +# Please follow this procedure in future cql3d modifications. + +SHELL = /bin/sh +NAME = xcql3d_mpi.gfortran64 +COMPILER= mpif77 +BUILDER= gfortran +INCLUDES = advnce.h comm.h frcomm.h frname.h frname_decl.h name.h \ + param.h trans.h wpadvnc.h name_decl.h +SOURCES = $(wildcard ADC/!(*mpitmp*).f) +SOURCES += a_cqlp.f abchief.f achief1.f achiefn.f aclear.f ainalloc.f \ + aindflt.f aindflt1.f aindfpa.f aingeom.f ainitial.f \ + ainpla.f ainplt.f ainpltpa.f \ + ainsetpa.f ainsetva.f ainspec.f ainvnorm.f ampfar.f \ + aminmx.f bavdens.f bavgmax.f baviorbt.f \ + bcast.f bsl.f bsu.f cfpcoefc.f cfpcoefn.f \ + cfpcoefr.f cfpgamma.f cfpleg.f \ + cfpmodbe.f cfpsymt.f coefefad.f coefefld.f coefegad.f \ + coeffpad.f coefload.f coefmidt.f coefmidv.f \ + coefrfad.f coefstup.f coefsyad.f coefwti.f \ + coefwtj.f diag.f diagcfac.f \ + diagdens.f diagdenz.f diagentr.f diagescl.f diaggnde.f \ + diaggnde2.f diagimpd.f diagscal.f diagwrng.f diagxswt.f \ + diagxswx.f dsk_gr.f dskout.f efield.f eflditer.f eqalloc.f \ + eqcoord.f eqelpse.f eqflxavg.f eqfn.f \ + eqfndpsi.f eqfninv.f eqfpsi.f \ + eqindflt.f eqinitl.f eqjac.f eqonovrp.f \ + eqorbit.f eqrhopsi.f eqrhs.f \ + eqtopeol.f equilib.f eqvolpsi.f eqwrng.f \ + esefld.f exlin.f exsweep.f exsweept.f \ + exsweepx.f finit.f firstdrv.f fle.f flxfn.f \ + freya.f freyasou.f frhexdrv.f \ + frinitl.f frinitz.f frnbdep2.f frnfreya.f \ + frplteq.f frset.f frsmooth.f frsplft.f \ + frstup.f frsubs.f frsuppor.f frwrong.f hpalloc0.f \ + ilut.f impavnc0.f impchk.f impnorm.f it3dalloc.f \ + lookup.f losscone.f lossegy.f lossorbm.f \ + losstor.f micfrplt.f micgetr.f \ + micgmbnd.f micgnbnd.f micxinil.f micxinim.f micxinit.f \ + micxiniz.f netcdfrf.f netcdfrw2.f ntdstore.f \ + ntloop.f pack21.f pltcont.f pltcycl.f pltdf.f \ + pltdnz.f pltelec.f pltendn.f \ + pltends.f pltfluxs.f pltfofvv.f pltfvsv.f pltinit.f pltlosc.f \ + pltmain.f pltpower.f pltprppr.f \ + pltrstv.f pltrun.f pltstrml.f pltvec.f pltvectr.f \ + pltvflux.f profaxis.f profiles.f prppr.f \ + prpprctr.f psif.f psifp.f psifppy.f psifpy.f psiinv.f \ + r8lsode.f r8subs.f rdc_multi.f rdc_bplt.f restcon.f resthks.f \ + restvty.f rf.f sigalloc.f siggy.f sigmax.f sigsetup.f \ + sigv5d.f sigfn.f sigie.f sigmaxwl.f sigv.f \ + soucrit.f sounorm.f soup.f soup0.f souplt.f \ + sourc0.f sourcee.f sourcef.f sourceko.f sourcpwr.f \ + synchrad.f tdbootst.f tdboothi.f tdchief.f tddiag.f \ + tddsig.f tdeqdsk.f tdfinterp.f tdinitl.f tdinlegw.f \ + tdinterp.f tdnflxs.f tdnpa.f tdnpadiag.f tdnpa0.f \ + tdnpacxcs.f tdnpalam.f tdnpabscs.f tdoutput.f \ + tdplteq.f tdpltjop.f tdpltmne.f tdpro.f \ + tdreadf.f tdrmshst.f tdsetnpa.f tdsetsxr.f \ + tdstin.f tdsxr.f tdsxr0.f tdsxray.f tdsxrplt.f \ + tdtloop.f tdtoarad.f \ + tdtoaray.f tdtraloc.f tdtransp.f tdtranspn.f tdtravct.f \ + tdtrchk.f tdtrchkd.f tdtrcon.f \ + tdtrdfus.f tdtrfcop.f tdtrflg.f tdtrflx.f \ + tdtrmuy.f tdtrrsou.f tdtrrtov.f \ + tdtrrtov2.f tdtrsavf.f tdtrsym.f tdtrvint.f \ + tdtrvsou.f tdtrvtor.f tdtrvtor2.f \ + tdtrvtor3.f tdtrwtl.f tdtry.f tdtscinp.f tdtscout.f \ + tdwrng.f tdwritef.f tdxin13d.f \ + tdxin23d.f tdxin33d.f tdxinitl.f urfalloc.f \ + urfavg.f urfb0.f urfbes.f \ + urfbplt.f urfchief.f urfdamp0.f urfdamp1.f \ + urfdamp2.f urfdampa.f urfdout.f urfedge.f \ + urffflx.f urfindfl.f urfinitl.f urfmidv.f \ + urfpack.f urfpackm.f urfrays.f urfread.f \ + urfread_.f urfsetup.f urfwrite.f urfwrite_.f \ + urfwrong.f urfwr0.f urfwr0c.f \ + vlf.f vlfalloc.f vlfbplt.f vlfsetup.f vlh.f vlhbplt.f vlhd.f \ + wpalloc.f wparsou.f wpavg.f wpbdry.f \ + wpcheck.f wpchgdy.f wpcthta.f wpelecf.f \ + wpinitl.f wploweq.f wpsavf.f wptrafx.f wptramu.f \ + wptrmuy.f wpvptb.f wpwrng.f wpmshchk.f \ + zblock.f zcunix.f zfreya.f znonsym.f +PYTHON = python3 +OBJECTS = $(SOURCES:.f=.o) +LOCATION = -L$(NETCDF_DIR)/lib -L$(PGPLOT_DIR)/lib -L$(PWD)/ADC +LIBS= -lX11 -lnetcdff -lnetcdf -lpgplot -lADC +INCLUDE=/usr/include +#DEBUG = -g -Mbounds +DEBUG = -g +OPTIMIZE = -O +LISTING = -Mlist +CSPECIAL = -Wno-padded +FFLAGS = $(CSPECIAL) $(OPTIMIZE) +#SPECIAL = -byteswapio -Mbackslash +LDSPECIAL = -Wl,-noinhibit-exec,-rpath,$(PGPLOT_DIR)/lib +#COMPILE = mpif77 -c $(CSPECIAL) -c $(DEBUG) -I $(INCLUDE) # or use $(OPTIMIZE) +COMPILE = $(COMPILER) -c $(FFLAGS) -I $(INCLUDE) # or use $(DEBUG) +LOAD = $(COMPILER) -o $(NAME) $(LDSPECIAL) # Remove -m for optimize +PROTECT = chmod 755 +DELETE = rm -f + +# The following gives suffixes to be used in checking for suffix rules. +# Written without dependencies, it may be useful to turn of such checking? +.SUFFIXES: + +$(NAME): $(OBJECTS) + $(LOAD) $(OBJECTS) $(LOCATION) $(LIBS) + $(PROTECT) $(NAME) + +# Following use of pattern matching works; +# it is based on statements in Sect. 10.7 of gmake Manual +# by Stallman and McGrath. +# However, the two rules following this seem clearer to me (BobH). +#%.o: %.f $(INCLUDES) +# $(COMPILE) $< -o $@ + +include deps #Introduced by John Wright (JCW) +$(SOURCES): $(INCLUDES) + +# see Sect. 4.10.1, Static Pattern Rules. +#$(OBJECTS): %.o: %.f $(INCLUDES) +%.o: %.f # $(INCLUDES) +#Changes introduced by JCW so obtain saved copy of mpi modified sources +# mpi/doparallel.py $< mpitmp.f mpi/mpins_par.f +# $(COMPILE) mpitmp.f -o $@ + export COMPILER FFLAGS + make -f ADC/Makefile + $(PYTHON) mpi/doparallel.py $< $*_mpitmp.f mpi/mpins_par.f + $(COMPILE) $*_mpitmp.f -o $@ + +rebuild: + $(COMPILE) $(SOURCES) + $(LOAD) $(OBJECTS) $(LOCATION) $(LIBS) + +clean: + make -f ADC/Makefile clean + $(DELETE) $(OBJECTS) *.lst *_mpitmp.f + diff --git a/mpi/doparallel.py b/mpi/doparallel.py index 04a7693..e6ca9a4 100644 --- a/mpi/doparallel.py +++ b/mpi/doparallel.py @@ -9,7 +9,7 @@ ########################################################################### -import string, sys +import sys, string inf = [] ## input lines out = [] ## output lines @@ -20,7 +20,7 @@ pat = [] ## pattern for cleaning buf = [] ## buffer for continuation lines -str = '' ## joined line +j_line = '' ## joined line ## clean single line def cleans(s): @@ -46,7 +46,7 @@ def cleans(s): elif state==2: if c=="'": state = 0 - if string.strip(s)=='': + if s.strip=='': return '' if len(r)<6: r += ' '*6 @@ -55,7 +55,7 @@ def cleans(s): ## saving buffer data into out list (c - prefix) def printbuf(c=''): global buf, out - if c!='' and string.strip(buf[0][:6])!='': + if c!='' and buf[0][:6].strip()!='': ## it is labeled line and it should be commented, so we use new line ## with given label and CONTINUE instruction out.append(buf[0][:6]+'CONTINUE !***') @@ -69,7 +69,7 @@ def printbuf(c=''): def procbuf(): if buf==[]: return - s = string.join(string.split(str),'') + s = j_line.strip() for p in pat: if p in s: printbuf('CMPI :::') @@ -78,11 +78,11 @@ def procbuf(): ## main cleaning module def clean(): - global buf, str, pat + global buf, j_line, pat ## reading patterns for cleaning P = open('mpi/patterns.mpi','r').readlines() for p in P: - pat.append(string.strip(p)) + pat.append(p.strip()) ## start processing buf = [] for s in inf: @@ -91,12 +91,12 @@ def clean(): buf.append(s) elif ss[5]!=' ': ## continuation line buf.append(s) - str += ss[6:] + j_line += ss[6:] continue else: ## normal line procbuf() buf = [s] - str = ss + j_line = ss procbuf() ########################################################################### @@ -111,7 +111,7 @@ def start(s, D): return False if len(s)>5 and s[5]!=' ': return False - c = string.strip(s[6:]) + c = s[6:].strip() if len(c)==0: return False pos = 0 @@ -121,7 +121,7 @@ def start(s, D): pos += 1 if pos==len(c): break - if string.lower(com) in D: + if com.lower() in D: return False return True @@ -136,7 +136,7 @@ def block(): ## checking that current file is in the list of blocked files check = False for l in open('mpi/block.mpi','r').readlines(): - if string.strip(l)==sys.argv[1]: + if l.strip()==sys.argv[1]: check = True break if not check: ## current file shouldn't be blocked @@ -144,7 +144,7 @@ def block(): ## reading delcaration instructions for detecting inserting point D = open('mpi/dec.mpi','r').readlines() for i in range(len(D)): - D[i] = string.strip(D[i]) + D[i] = D[i].strip() ## searching insertion point and inserting blocking instruction state = 0 for i in range(len(out)): @@ -165,7 +165,7 @@ def dompi(): L = open(sys.argv[3],'r').readlines() for l in L: if l[0]=='C': - key = string.strip(l) + key = l.strip() ins[key] = [] else: ins[key].append(l[:-1]) @@ -176,8 +176,8 @@ def dompi(): inf.append('CMPI :::'+l) mode = '' continue - ll = string.strip(l) - if l[0:4]=='CMPI' and ins.has_key(ll): + ll = l.strip() + if l[0:4]=='CMPI' and ll in ins: inf.append('CMPI >>>') for i in ins[ll]: inf.append(i) @@ -195,10 +195,14 @@ def main(): global inf, out ## reading input file into inf buffer + if "mpilib.f" in sys.argv[1]: + # Don't process this source file as its contents will be parsed + # into others + return f = open(sys.argv[1],'r') s = f.readline() while s!='': - inf.append(string.rstrip(s)) + inf.append(s.rstrip()) s = f.readline() f.close() @@ -210,7 +214,7 @@ def main(): ## writing result f = open(sys.argv[2],'w') for s in inf: - print >> f, s + f.write(s + "\r\n") f.close() ## do all diff --git a/pltinit.f b/pltinit.f index cfb16ad..3c14311 100644 --- a/pltinit.f +++ b/pltinit.f @@ -26,7 +26,7 @@ subroutine pltinit c ier=pgbeg(0,'plot.ps/VPS',1,1) c Remember, pgopen should be integer write(t_,1000) mnemonic(1:length_char(mnemonic)) - 1000 format(a,".ps/VCPS") !YuP: was /VPS (vertical black&white) + 1000 format(a,".ps/PS") !YuP: was /VPS (vertical black&white) ! Use /VCPS for vertical Color pages ier=PGOPEN(t_) CALL PGSCI(1) @@ -58,14 +58,14 @@ subroutine pltinit !'CURSOR' * - does the current device have a graphics cursor? ! ('YES' or 'NO'). ! Two other arg. are outputs: - CALL PGQINF('TYPE', PG_VAL, PG_L) - WRITE (*,*) 'PGPLOT device type: ', PG_VAL(1:PG_L) - CALL PGQINF('DEVICE', PG_VAL, PG_L) - WRITE (*,*) 'PGPLOT device: ', PG_VAL(1:PG_L) - CALL PGQINF('USER', PG_VAL, PG_L) - WRITE (*,*) 'PGPLOT user: ', PG_VAL(1:PG_L) - CALL PGQINF('NOW', PG_VAL, PG_L) - WRITE (*,*) 'PGPLOT time now: ', PG_VAL(1:PG_L) + ! CALL PGQINF('TYPE', PG_VAL, PG_L) + ! WRITE (*,*) 'PGPLOT device type: ', PG_VAL(1:PG_L) + ! CALL PGQINF('DEVICE', PG_VAL, PG_L) + ! WRITE (*,*) 'PGPLOT device: ', PG_VAL(1:PG_L) + ! CALL PGQINF('USER', PG_VAL, PG_L) + ! WRITE (*,*) 'PGPLOT user: ', PG_VAL(1:PG_L) + ! CALL PGQINF('NOW', PG_VAL, PG_L) + ! WRITE (*,*) 'PGPLOT time now: ', PG_VAL(1:PG_L) !Inquire color index range: CALL PGQCIR(PG_C1, PG_C2)