! $Id: aerosol_init.f90 21410 2013-12-25 15:44:43Z Nbabkovskaia ! ! This module provide a way for users to specify custom initial ! conditions. ! ! The module provides a set of standard hooks into the Pencil Code ! and currently allows the following customizations: ! ! Description | Relevant function call ! ------------------------------------------------------------------------ ! Initial condition registration | register_initial_condition ! (pre parameter read) | ! Initial condition initialization | initialize_initial_condition ! (post parameter read) | ! | ! Initial condition for momentum | initial_condition_uu ! Initial condition for density | initial_condition_lnrho ! Initial condition for entropy | initial_condition_ss ! Initial condition for magnetic potential | initial_condition_aa ! ! And a similar subroutine for each module with an "init_XXX" call. ! The subroutines are organized IN THE SAME ORDER THAT THEY ARE CALLED. ! First uu, then lnrho, then ss, then aa, and so on. ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! CPARAM logical, parameter :: linitial_condition = .true. ! !*************************************************************** module InitialCondition ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages use EquationOfState ! implicit none ! include '../initial_condition.h' ! real :: init_ux=impossible,init_uy=impossible,init_uz=impossible integer :: spot_number=10 integer :: index_H2O=2 integer :: index_N2=3, i_point=1 integer :: Ndata=10, Nadd_points=0 real :: dYw=1.,dYw1=1.,dYw2=1., init_water1=0., init_water2=0. real :: init_x1=0.,init_x2=0.,init_TT1=0., init_TT2=0. real, dimension(nchemspec) :: init_Yk_1, init_Yk_2 real :: X_wind=impossible, spot_size=10. real :: AA=0.66e-4, d0=2.4e-6 , BB0=1.5*1e-16, rhow_coeff=1. real :: dsize_min=0., dsize_max=0., r0=0., r02=0., Period=2., delta real, dimension(ndustspec) :: dsize, dsize0 real, dimension(20000) :: Ntot_data ! real, dimension(mx,ndustspec) :: init_distr_loc logical :: lreinit_water=.false.,lwet_spots=.false. logical :: linit_temperature=.false., lcurved_xz=.false. logical :: ltanh_prof_xy=.false.,ltanh_prof_xz=.false. logical :: llog_distribution=.true., lcurved_xy=.false. logical :: lACTOS=.false.,lACTOS_read=.true., lACTOS_write=.true., lsinhron=.false. logical :: ladd_points=.false., lrho_const=.false., lregriding=.false. logical :: lLES=.false., lP_aver=.false. ! namelist /initial_condition_pars/ & init_ux, init_uy,init_uz,init_x1,init_x2, init_water1, init_water2, & lreinit_water, dYw,dYw1, dYw2, X_wind, spot_number, spot_size, lwet_spots, & linit_temperature, init_TT1, init_TT2, dsize_min, dsize_max, r0, r02, d0, lcurved_xz, lcurved_xy, & ltanh_prof_xz,ltanh_prof_xy, Period, BB0, index_N2, index_H2O, lACTOS, lACTOS_read, lACTOS_write, & i_point,Ndata, lsinhron, delta, Nadd_points, ladd_points, lrho_const, lregriding, lLES, lP_aver, rhow_coeff ! contains !*********************************************************************** subroutine register_initial_condition() ! ! Register variables associated with this module; likely none. ! ! 07-may-09/wlad: coded ! if (lroot) call svn_id( & "$Id$") ! endsubroutine register_initial_condition !*********************************************************************** subroutine initialize_initial_condition(f) ! ! Initialize any module variables which are parameter dependent. ! ! 07-may-09/wlad: coded ! real, dimension (mx,my,mz,mfarray) :: f ! call keep_compiler_quiet(f) ! endsubroutine initialize_initial_condition !*********************************************************************** subroutine initial_condition_uu(f) ! ! Initialize the velocity field. ! ! 07-may-09/wlad: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f integer :: i,j real :: del=10.,bs ! if ((init_ux /=impossible) .and. (nygrid>1)) then do i=1,my f(:,i,:,iux)=cos(Period*PI*y(i)/Lxyz(2))*init_ux enddo endif ! if ((init_uz /=impossible) .and. (nzgrid>1)) then ! bs=9.81e2*(293.-290.)/293. ! do i=1,mx ! f(i,:,:,iuz)=-sqrt(Lxyz(1)*bs) & ! *((exp(2.*x(i)/Lxyz(1))+exp(-2.*x(i)/Lxyz(1)))/2.)**(-2) ! enddo ! endif ! if ((init_uy /=impossible) .and. (X_wind /= impossible)) then do j=1,mz f(:,:,j,iuy)= init_uy*(z(j)-xyz0(3))/(Lxyz(3)-xyz0(3)) ! +(init_uy+0.)*0.5+((init_uy-0.)*0.5) & ! *(exp((x(j)+X_wind)/del)-exp(-(x(j)+X_wind)/del)) & ! /(exp((x(j)+X_wind)/del)+exp(-(x(j)+X_wind)/del)) enddo endif ! if ((init_ux /=impossible) .and. (X_wind /= impossible)) then do j=1,mx f(:,:,j,iux)=init_ux*(z(j)-xyz0(3))/(Lxyz(3)-xyz0(3)) ! ! +(init_uz+0.)*0.5+((init_uz-0.)*0.5) & ! *(exp((x(j)+X_wind)/del)-exp(-(x(j)+X_wind)/del)) & ! /(exp((x(j)+X_wind)/del)+exp(-(x(j)+X_wind)/del)) enddo endif ! if ((init_uy /=impossible) .and. (X_wind == impossible)) then f(:,:,:,iuy)=f(:,:,:,iuy)+init_uy endif if ((init_uz /=impossible) .and. (X_wind == impossible)) then ! ! f(:,:,:,iuz)=f(:,:,:,iuz)+init_uz ! do i=l1,l2 if (x(i)>0) then f(i,:,:,iuz)=init_uz else f(i,:,:,iuz)=-init_uz endif enddo ! endif if ((init_ux /=impossible) .and. (X_wind == impossible)) then f(:,:,:,iux)=f(:,:,:,iux)+init_ux endif ! call keep_compiler_quiet(f) ! endsubroutine initial_condition_uu !*********************************************************************** subroutine initial_condition_lnrho(f) ! ! Initialize logarithmic density. init_lnrho will take care of ! converting it to linear density if you use ldensity_nolog. ! ! 07-may-09/wlad: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! ! SAMPLE IMPLEMENTATION ! call keep_compiler_quiet(f) ! endsubroutine initial_condition_lnrho !*********************************************************************** subroutine initial_condition_ss(f) ! ! Initialize entropy. ! ! 07-may-09/wlad: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! ! SAMPLE IMPLEMENTATION ! call keep_compiler_quiet(f) ! endsubroutine initial_condition_ss !*********************************************************************** subroutine initial_condition_chemistry(f) ! ! Initialize chemistry. ! ! 07-may-09/wlad: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f real, dimension (my,mz) :: init_water1,init_water2, & init_water1_tmp,init_water2_tmp !real, dimension (mx,my,mz), intent(inout) :: f real, dimension (ndustspec) :: lnds real :: ddsize, ddsize0, del, air_mass, PP, TT integer :: i, ii logical :: lstart1=.false., lstart2=.false. ! if (llog_distribution) then ddsize=(alog(dsize_max)-alog(dsize_min))/(max(ndustspec,2)-1) else ddsize=(dsize_max-dsize_min)/(max(ndustspec,2)-1) endif ! do i=0,(ndustspec-1) if (llog_distribution) then lnds(i+1)=alog(dsize_min)+i*ddsize dsize(i+1)=exp(lnds(i+1)) else lnds(i+1)=dsize_min+i*ddsize dsize(i+1)=lnds(i+1) endif enddo ! if (lACTOS) then call ACTOS_data(f) call reinitialization(f, air_mass, PP, TT) elseif (lLES) then call LES_data(f) else call air_field_local(f, air_mass, PP, TT) call reinitialization(f, air_mass, PP, TT) endif ! ! call reinitialization(f, air_mass, PP, TT) ! endsubroutine initial_condition_chemistry !*********************************************************************** subroutine initial_condition_uud(f) ! ! Initialize dust fluid velocity. ! ! 07-may-09/wlad: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine initial_condition_uud !*********************************************************************** subroutine initial_condition_nd(f) ! ! Initialize dust fluid density. ! ! 07-may-09/wlad: coded ! use General, only: polynomial_interpolation real, dimension (mx,my,mz,mfarray), intent(inout) :: f integer :: i_spline, k, i real, dimension (4) :: ary, arx real, dimension (3) :: x2, S, ddy real, dimension (ndustspec) :: nd_tmp real :: tmp2, ddsize real, dimension(ndustspec) :: lnds, dsize ! call keep_compiler_quiet(f) ! endsubroutine initial_condition_nd !*********************************************************************** subroutine initial_condition_uun(f) ! ! Initialize neutral fluid velocity. ! ! 07-may-09/wlad: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine initial_condition_uun !*********************************************************************** subroutine read_initial_condition_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=initial_condition_pars, IOSTAT=iostat) ! endsubroutine read_initial_condition_pars !*********************************************************************** subroutine write_initial_condition_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=initial_condition_pars) ! endsubroutine write_initial_condition_pars !*********************************************************************** subroutine air_field_local(f, air_mass, PP, TT) ! real, dimension (mx,my,mz,mvar+maux) :: f real, dimension (mx,my,mz) :: sum_Y, tmp ! logical :: emptyfile=.true. logical :: found_specie integer :: file_id=123, ind_glob, ind_chem character (len=80) :: ChemInpLine character (len=10) :: specie_string character (len=1) :: tmp_string integer :: i,j,k=1,index_YY, j1,j2,j3, iter real :: YY_k, air_mass real, intent(out) :: PP, TT ! (in dynes = 1atm) real, dimension(nchemspec) :: stor2 integer, dimension(nchemspec) :: stor1 ! integer :: StartInd,StopInd,StartInd_1,StopInd_1 integer :: iostat, i1,i2,i3 ! air_mass=0. StartInd_1=1; StopInd_1 =0 open(file_id,file="air.dat") ! if (lroot) print*, 'the following parameters and '//& 'species are found in air.dat (volume fraction fraction in %): ' ! dataloop: do ! read(file_id,'(80A)',IOSTAT=iostat) ChemInpLine if (iostat < 0) exit dataloop emptyFile=.false. StartInd_1=1; StopInd_1=0 StopInd_1=index(ChemInpLine,' ') specie_string=trim(ChemInpLine(1:StopInd_1-1)) tmp_string=trim(ChemInpLine(1:1)) ! if (tmp_string == '!' .or. tmp_string == ' ') then elseif (tmp_string == 'T') then StartInd=1; StopInd =0 ! StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 StartInd=verify(ChemInpLine(StopInd:),' ')+StopInd-1 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 ! read (unit=ChemInpLine(StartInd:StopInd),fmt='(E14.7)') TT if (lroot) print*, ' Temperature, K ', TT ! elseif (tmp_string == 'P') then ! StartInd=1; StopInd =0 ! StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 StartInd=verify(ChemInpLine(StopInd:),' ')+StopInd-1 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 ! read (unit=ChemInpLine(StartInd:StopInd),fmt='(E14.7)') PP if (lroot) print*, ' Pressure, Pa ', PP ! else ! call find_species_index(specie_string,ind_glob,ind_chem,found_specie) ! if (found_specie) then ! StartInd=1; StopInd =0 ! StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 StartInd=verify(ChemInpLine(StopInd:),' ')+StopInd-1 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') YY_k if (lroot) print*, ' volume fraction, %, ', YY_k, & species_constants(ind_chem,imass) ! if (species_constants(ind_chem,imass)>0.) then air_mass=air_mass+YY_k*0.01/species_constants(ind_chem,imass) endif ! if (StartInd==80) exit ! stor1(k)=ind_chem stor2(k)=YY_k k=k+1 endif ! endif enddo dataloop ! ! Stop if air.dat is empty ! if (emptyFile) call fatal_error("air_field", "I can only set existing fields") air_mass=1./air_mass ! do j=1,k-1 f(:,:,:,ichemspec(stor1(j)))=stor2(j)*0.01 enddo ! sum_Y=0. ! do j=1,nchemspec sum_Y=sum_Y+f(:,:,:,ichemspec(j)) enddo do j=1,nchemspec f(:,:,:,ichemspec(j))=f(:,:,:,ichemspec(j))/sum_Y enddo ! ! do j=1,nchemspec init_Yk_1(j)=f(l1,m1,n1,ichemspec(j)) init_Yk_2(j)=f(l1,m1,n1,ichemspec(j)) enddo ! if (mvar < 5) then call fatal_error("air_field", "I can only set existing fields") endif if (ltemperature_nolog) then f(:,:,:,iTT)=TT else f(:,:,:,ilnTT)=alog(TT)!+f(:,:,:,ilnTT) endif if (ldensity_nolog) then f(:,:,:,ilnrho)=(PP/(k_B_cgs/m_u_cgs)*& air_mass/TT)/unit_mass*unit_length**3 else tmp=(PP/(k_B_cgs/m_u_cgs)*& air_mass/TT)/unit_mass*unit_length**3 f(:,:,:,ilnrho)=alog(tmp) endif ! if (ltemperature_nolog) then f(:,:,:,iTT)=TT else f(:,:,:,ilnTT)=alog(TT)!+f(:,:,:,ilnTT) endif if (ldensity_nolog) then f(:,:,:,ilnrho)=(PP/(k_B_cgs/m_u_cgs)*& air_mass/TT)/unit_mass*unit_length**3 else tmp=(PP/(k_B_cgs/m_u_cgs)*& air_mass/TT)/unit_mass*unit_length**3 f(:,:,:,ilnrho)=alog(tmp) endif ! if (lroot) print*, 'local:Air temperature, K', TT if (lroot) print*, 'local:Air pressure, dyn', PP if (lroot) print*, 'local:Air density, g/cm^3:' if (lroot) print '(E10.3)', PP/(k_B_cgs/m_u_cgs)*air_mass/TT if (lroot) print*, 'local:Air mean weight, g/mol', air_mass if (lroot) print*, 'local:R', k_B_cgs/m_u_cgs ! close(file_id) !! endsubroutine air_field_local !*********************************************************************** subroutine ACTOS_data(f) ! real, dimension (mx,my,mz,mvar+maux) :: f real, dimension (mx,my,mz) :: sum_Y, tmp, air_mass real, dimension (20000) :: PP_data, rhow_data, TT_data, tmp_data real, dimension (20000) :: PP_data_add, rhow_data_add, TT_data_add real, dimension (20000) :: ttime2 real, dimension (mx,my,mz) :: ux_data, uy_data, uz_data real, dimension (1340) :: ttime real, dimension (1300,6) :: coeff_loc real, dimension (20000,7) :: coeff_loc2 real, dimension (20000,ndustspec) :: init_distr_loc, init_distr_tmp real, dimension (6) :: tmp3 real, dimension (ndustspec) :: tmp4 ! logical :: emptyfile=.true., lfind logical :: found_specie integer :: file_id=123, ind_glob, ind_chem,jj character (len=800) :: ChemInpLine integer :: i,j,k=1,index_YY, j1,j2,j3, iter, ll1, mm1, nn1 real :: TT=300., ddsize, tmp2, right, left, PP_aver=0., rho_aver double precision, dimension (mx,my,mz) :: tmp5 ! real, intent(out) :: PP ! (in dynes = 1atm) real, dimension(nchemspec) :: stor2, stor1 real, dimension(29) :: input_data, input_data2 ! real, dimension (7) :: ctmp integer :: StartInd,StopInd,StartInd_1,StopInd_1 integer :: iostat, i1,i2,i3, i1_left,i1_right logical :: lwrite_string=.false. ! if (lACTOS_write) then ! if (lsinhron) then ! open(143,file="coeff_part.dat") ! do i=1,1300 ! read(143,'(f15.6,f15.6,f15.6,f15.6,f15.6,f15.6,f15.6)'),tmp3,ttime(i) ! coeff_loc(i,:)=tmp3 ! enddo ! close(143) ! endif ! ! air_mass=0. StartInd_1=1; StopInd_1 =0 open(file_id,file="ACTOS_data.out") open(143,file="ACTOS_new.out") open(144,file="coeff_part_new.out") ! open(file_id,file="ACTOS_xyz_data.out") ! open(143,file="ACTOS_xyz_new.out") ! j=1 i=1 lwrite_string=.true. do while (j<=780000) read(file_id,'(80A)',IOSTAT=iostat) ChemInpLine StartInd=1; StopInd =0 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 ! if ((i>1) .and. (i_point>1)) then if (i1) .and.(i_point==1)) then lwrite_string=.true. endif !------------------------------------------ if (lwrite_string) then k=1 do while (k<30) if (StopInd==StartInd) then StartInd=StartInd+1 else if (k<29) then StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 else StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 endif read (unit=ChemInpLine(StartInd:StopInd-1),fmt='(f12.6)') tmp2 input_data(k)=tmp2 StartInd=StopInd k=k+1 endif enddo ! !***************************************** if (lsinhron) then ttime2(j)=input_data(23) lfind=.false. do jj=1,1300 if (abs(ttime(jj)-input_data(23))<.05) then lfind=.true. else lfind=.false. endif if ((ttime2(j) == ttime2(j-1)) .or. (j==1)) lfind=.false. if (lfind) then write(143,'(29f15.6)') input_data print*, jj, input_data(23), ttime(jj), lfind, abs(ttime(jj)-input_data(23)) endif enddo else ! lfind=.false. !! this is full data if (input_data(23)>52000.) then !! this is data for calculations !! if ((input_data(23)>54200.) .and. (input_data(23)<54240.)) then write(143,'(29f15.6)') input_data endif endif !************************************** i=i+1 lwrite_string=.false. endif !-------------------------------------- j=j+1 enddo ! close(file_id) close(143) close(144) ! endif ! ! emptyFile=.false. ! ! if (lACTOS_read) then ! open(143,file="ACTOS_new.out") do i=1,Ndata read(143,'(29f15.6)') input_data TT_data(i)=input_data(10)+272.15 PP_data(i)=input_data(7)*1e3 !dyn rhow_data(i)=input_data(16)*1e-6 !g/cm3 ! uvert(i)=input_data(26)*1e2 ! ux_data(i)=input_data(25)*1e2 ! uy_data(i)=input_data(26)*1e2 ! uz_data(i)=input_data(27)*1e2 ! Ntot_data(i)=input_data(11) enddo close(143) ! if (lregriding) then open(143,file="ux.dat") do k=1,mz do j=1,my do i=1,mx read(143,'(e9.5)') tmp2 ux_data(i,j,k)=tmp2 f(i,j,k,iux)=ux_data(i,j,k) enddo enddo enddo close(143) endif ! ! open(143,file="ACTOS_xyz_new.out") ! do i=1,Ndata ! read(143,'(29f15.6)') input_data2 ! ux_data(i)=uvert(i)/cos(input_data2(19))/cos(input_data2(20)) ! uy_data(i)=uvert(i)/cos(input_data2(19))/sin(input_data2(20)) ! uz_data(i)=uvert(i)/sin(input_data2(19)) ! enddo ! close(143) ! PP_aver=0. do i=1,Ndata PP_aver=PP_aver+PP_data(i) enddo PP_data=PP_aver/Ndata print*,'PP_aver=',PP_aver/Ndata ! ! print*,TT_data(int(0.05*nygrid)),TT_data(nygrid-int(0.05*nygrid)) ! if (ladd_points) then k=1 do i=1,Ndata TT_data_add(k)=TT_data(i) PP_data_add(k)=PP_data(i) rhow_data_add(k)=rhow_data(i) k=k+1 do j=1,Nadd_points TT_data_add(k)=TT_data(i)+(TT_data(i+1)-TT_data(i))*j/(Nadd_points+1) PP_data_add(k)=PP_data(i)+(PP_data(i+1)-PP_data(i))*j/(Nadd_points+1) rhow_data_add(k)=rhow_data(i)+(rhow_data(i+1)-rhow_data(i))*j/(Nadd_points+1) k=k+1 ! print*,'k=',k enddo enddo endif mm1=anint((y(l1)-xyz0(2))/dy) do i=m1,m2 if (ladd_points) then f(:,i,:,ilnTT)=alog(TT_data_add(mm1+i-3)) f(:,i,:,ichemspec(index_H2O))=rhow_data_add(mm1+i-3)/1e-2/10. !g/cm3 else f(:,i,:,ilnTT)=alog(TT_data(mm1+i-3)) f(:,i,:,ichemspec(index_H2O))=rhow_data(mm1+i-3)/1e-2/10. !g/cm3 endif enddo ! if (nygrid>1) then ! mm1=anint((y(m1)-xyz0(2))/dy) ! ! if (init_uy /= impossible) then ! f(:,:,:,iuy)=init_uy ! else ! do i=m1,m2 ! f(:,i,:,iuy)=uy_data(mm1+i-3) ! enddo ! endif ! ! endif ! if (nzgrid>1) then ! nn1=anint((z(m1)-xyz0(3))/dz) ! ! if (init_uz /= impossible) then ! f(:,:,:,iuz)=init_uz ! else ! do i=n1,n2 ! f(:,:,i,iuz)=uz_data(nn1+i-3) ! enddo ! endif ! ! endif ! f(:,:,:,ichemspec(index_N2))=0.7 f(:,:,:,ichemspec(1))=1.-f(:,:,:,ichemspec(index_N2))-f(:,:,:,ichemspec(index_H2O)) ! Stop if air.dat is empty ! ! if (emptyFile) call fatal_error("ACTOS data", "I can only set existing fields") ! air_mass=1./air_mass ! ! if (lrho_const) then sum_Y=0. do k=1,nchemspec sum_Y=sum_Y + f(:,:,:,ichemspec(k))/species_constants(k,imass) enddo air_mass=1./sum_Y ! do i=1,Ndata tmp_data(i)=PP_data(i)/TT_data(i) enddo ! rho_aver=sum(tmp_data)/Ndata*sum(air_mass)/mx/my/mz rho_aver=rho_aver/(k_B_cgs/m_u_cgs)/unit_mass*unit_length**3 endif do iter=1,4 ! sum_Y=0. do k=1,nchemspec sum_Y=sum_Y + f(:,:,:,ichemspec(k))/species_constants(k,imass) enddo air_mass=1./sum_Y ! do i=m1,m2 if (ladd_points) then tmp5(:,i,:)=dlog(PP_data_add(mm1+i-3)/(k_B_cgs/m_u_cgs)*air_mass(:,i,:) & /exp(f(:,i,:,ilnTT))/unit_mass*unit_length**3) else tmp5(:,i,:)=dlog(PP_data(mm1+i-3)/(k_B_cgs/m_u_cgs)*air_mass(:,i,:) & /exp(f(:,i,:,ilnTT))/unit_mass*unit_length**3) endif if (lrho_const) then f(:,i,:,ilnrho)=alog(rho_aver) else f(:,i,:,ilnrho)=tmp5(:,i,:) endif enddo ! if (iter<4) then do i=m1,m2 if (ladd_points) then f(:,i,:,ichemspec(index_H2O))=rhow_data_add(mm1+i-3)/exp(f(:,i,:,ilnrho))*1. else f(:,i,:,ichemspec(index_H2O))=rhow_data(mm1+i-3)/exp(f(:,i,:,ilnrho))*1. endif enddo f(:,:,:,ichemspec(1))=1.-f(:,:,:,ichemspec(index_N2))-f(:,:,:,ichemspec(index_H2O)) endif ! enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! particles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! open(143,file="coeff_part_new.out") ! do i=1,Ndata ! read(143,'(f15.6,f15.6,f15.6,f15.6,f15.6,f15.6,f15.6)'),ctmp ! coeff_loc2(i,:)=ctmp ! enddo ! close(143) ! ! do i=1,Ndata ! do k=1,ndustspec ! tmp2=dsize(k)*1e4*2. ! if (coeff_loc2(i,3)/=0.) then ! init_distr_loc(i,k) = (coeff_loc2(i,1) & ! *exp(-0.5*( (tmp2-coeff_loc2(i,2)) /coeff_loc2(i,3) )**2) & ! +coeff_loc2(i,4)+coeff_loc2(i,5)*tmp2+coeff_loc2(i,6)*tmp2**2) ! else ! init_distr_loc(i,k)=0. ! endif ! if (init_distr_loc(i,k)<0) init_distr_loc(i,k)=0. ! enddo ! enddo ! ! do k=1,ndustspec ! do i=1,Ndata ! if (init_distr_loc(i,k)==0.) then ! i1=0; right=0.; i1_right=0 ! do while ((right == 0.) .and. (i1_right0)) ! left=init_distr_loc(i-i1,k) ! i1_left=i-i1 ! i1=i1+1 ! enddo ! ! if (i1_left==1) left=right ! if (i1_right==Ndata) right=left ! ! if ((right==0.) .and. (left==0.)) then ! init_distr_tmp(i,k)=0. ! else ! init_distr_tmp(i,k)=left+(right-left) & ! *(coeff_loc2(i,7)-coeff_loc2(i1_left,7))/(coeff_loc2(i1_right,7)-coeff_loc2(i1_left,7)) ! endif ! endif ! enddo ! enddo ! ! do k=1,ndustspec ! do i=1,Ndata ! if (init_distr_loc(i,k)==0.) init_distr_loc(i,k)=init_distr_tmp(i,k) ! if ((2.*dsize(k)*1e4<.5) .or. (2.*dsize(k)*1e4>40.)) init_distr_loc(i,k)=0. ! enddo ! enddo ! ! open(144,file="part_new.out") ! do i=1,Ndata ! tmp4=init_distr_loc(i,:) ! write(144,'(29f15.6)'),coeff_loc2(i,7),tmp4 ! enddo ! close(144) ! ! do i=l1,l2 ! do k=1,ndustspec ! f(i,:,:,ind(k)) = (Ntot_data(ll1+i-3)/(2.*pi)**0.5/alog(delta) & ! (init_distr_loc(ll1+i-3,k) + Ntot_data(ll1+i-3)/(2.*pi)**0.5/alog(delta) & ! * exp(-(alog(2.*dsize(k))-alog(2.*r0))**2/(2.*(alog(delta))**2))) & ! /exp(f(i,:,:,ilnrho))/dsize(k) ! enddo ! enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! if (lroot) print*, 'local:Air temperature, K', maxval(exp(f(l1:l2,m1:m2,n1:n2,ilnTT))), & minval(exp(f(l1:l2,m1:m2,n1:n2,ilnTT))) if (lroot) print*, 'local:Air pressure, dyn', maxval(PP_data), minval(PP_data) if (lroot) print*, 'local:Air density, g/cm^3:' if (lroot) print '(E10.3)', maxval(exp(f(:,:,:,ilnrho))), minval(exp(f(:,:,:,ilnrho))) if (lroot) print*, 'local:Air mean weight, g/mol', maxval(air_mass),minval(air_mass) if (lroot) print*, 'local:R', k_B_cgs/m_u_cgs ! endif ! endsubroutine ACTOS_data !*********************************************************************** subroutine LES_data(f) ! real, dimension (mx,my,mz,mvar+maux) :: f real, dimension (mx,my,mz) :: sum_Y, tmp, air_mass real, dimension (2000) :: PP_data, rhow_data, TT_data, tmp_data, w_data ! real, dimension (20000) :: ttime2 real, dimension (mx,my,mz) :: ux_data, uy_data, uz_data ! real, dimension (1340) :: ttime ! real, dimension (1300,6) :: coeff_loc ! real, dimension (20000,7) :: coeff_loc2 real, dimension (20000,ndustspec) :: init_distr_loc, init_distr_tmp real, dimension (6) :: tmp3 real, dimension (ndustspec) :: tmp4 ! logical :: emptyfile=.true., lfind logical :: found_specie integer :: file_id=123, ind_glob, ind_chem,jj character (len=800) :: ChemInpLine integer :: i,j,k=1,index_YY, iter, nn1, ii, io_code real :: TT=300., ddsize, tmp2, right, left, PP_aver=0., rho_aver double precision, dimension (mx,my,mz) :: tmp5 ! real, intent(out) :: PP ! (in dynes = 1atm) real, dimension(nchemspec) :: stor2, stor1 real, dimension(2) :: input_data ! real, dimension (7) :: ctmp integer :: StartInd,StopInd,StartInd_1,StopInd_1 integer :: iostat, i1,i2,i3, i1_left,i1_right logical :: lwrite_string=.false. ! open(143,file="T_init.dat") do i=1,Ndata read(143,*,iostat=io_code) (input_data(ii),ii=1,2) TT_data(i)=input_data(2) ! print*,i,' ',TT_data(i) enddo close(143) ! open(143,file="pre_init.dat") do i=1,Ndata read(143,*,iostat=io_code) (input_data(ii),ii=1,2) PP_data(i)=input_data(2)*10. !dyn enddo close(143) ! open(143,file="qv_init.dat") do i=1,Ndata read(143,*,iostat=io_code) (input_data(ii),ii=1,2) rhow_data(i)=input_data(2) !dyn ! print*,i,' ',rhow_data(i) enddo close(143) open(143,file="wind.dat") do i=1,Ndata read(143,*,iostat=io_code) (input_data(ii),ii=1,2) w_data(i)=input_data(2) !dyn ! ! print*,i,' ',w_data(i) enddo close(143) if (lP_aver) then PP_aver=0. do i=1,Ndata PP_aver=PP_aver+PP_data(i) enddo PP_data=PP_aver/Ndata endif ! ! print*,'PP_aver=',PP_aver/Ndata ! ! print*,TT_data(int(0.05*nygrid)),TT_data(nygrid-int(0.05*nygrid)) ! nn1=anint((z(n1)-xyz0(3))/dz) do i=n1,n2 f(:,:,i,ilnTT)=alog(TT_data(nn1+i-3)) f(:,:,i,ichemspec(index_H2O))=rhow_data(nn1+i-3)*rhow_coeff enddo ! f(:,:,:,ichemspec(index_N2))=0.7 f(:,:,:,ichemspec(1))=1.-f(:,:,:,ichemspec(index_N2))-f(:,:,:,ichemspec(index_H2O)) ! print*,maxval(f(:,:,:,ichemspec(1))), minval(f(:,:,:,ichemspec(1))) ! Stop if air.dat is empty ! ! if (emptyFile) call fatal_error("ACTOS data", "I can only set existing fields") ! if (lrho_const) then sum_Y=0. do k=1,nchemspec sum_Y=sum_Y + f(:,:,:,ichemspec(k))/species_constants(k,imass) enddo air_mass=1./sum_Y ! do i=1,Ndata tmp_data(i)=PP_data(i)/TT_data(i) enddo ! rho_aver=sum(tmp_data)/Ndata*sum(air_mass)/mx/my/mz rho_aver=rho_aver/(k_B_cgs/m_u_cgs)/unit_mass*unit_length**3 endif do iter=1,4 ! sum_Y=0. do k=1,nchemspec sum_Y=sum_Y + f(:,:,:,ichemspec(k))/species_constants(k,imass) ! print*,k,' ',species_constants(k,imass) enddo air_mass=1./sum_Y ! do i=n1,n2 tmp5(:,:,i)=dlog(PP_data(nn1+i-3)/(k_B_cgs/m_u_cgs)*air_mass(:,:,i) & /exp(f(:,:,i,ilnTT))/unit_mass*unit_length**3) if (lrho_const) then f(:,:,i,ilnrho)=alog(rho_aver) else f(:,:,i,ilnrho)=tmp5(:,:,i) endif f(:,:,i,iux)=w_data(nn1+i-3)*100.*sqrt(3.)/2. f(:,:,i,iuy)=w_data(nn1+i-3)*100.*1./2. ! print*,w_data(nn1+i-3),i, nn1+i-3 enddo ! ! if (iter<4) then ! do i=n1,n2 ! f(:,:,i,ichemspec(index_H2O))=rhow_data(nn1+i-3)/exp(f(:,:,i,ilnrho))*1. ! enddo ! f(:,:,:,ichemspec(1))=1.-f(:,:,:,ichemspec(index_N2))-f(:,:,:,ichemspec(index_H2O)) ! endif ! enddo ! if (lroot) print*, 'local:Air temperature, K', maxval(exp(f(l1:l2,m1:m2,n1:n2,ilnTT))), & minval(exp(f(l1:l2,m1:m2,n1:n2,ilnTT))) if (lroot) print*, 'local:Air pressure, dyn', maxval(PP_data), minval(PP_data) if (lroot) print*, 'local:Air density, g/cm^3:' if (lroot) print '(E10.3)', maxval(exp(f(:,:,:,ilnrho))), minval(exp(f(:,:,:,ilnrho))) if (lroot) print*, 'local:Air mean weight, g/mol', maxval(air_mass),minval(air_mass) if (lroot) print*, 'local:R', k_B_cgs/m_u_cgs ! ! print*,'2',maxval(f(:,:,:,ichemspec(1))), minval(f(:,:,:,ichemspec(1))) endsubroutine LES_data !*********************************************************************** subroutine reinitialization(f, air_mass, PP, TT) ! real, dimension (mx,my,mz,mvar+maux) :: f real, dimension (mx,my,mz) :: sum_Y, air_mass_ar, tmp real, dimension (mx,my,mz) :: init_water1_,init_water2_ real, dimension (my,mz) :: init_water1_min,init_water2_max real , dimension (my) :: init_x1_ary, init_x2_ary, del_ary, del_ar1y, del_ar2y real , dimension (mz) :: init_x1_arz, init_x2_arz, del_arz, del_ar1z, del_ar2z ! integer :: i,j,k, j1,j2,j3, iter real :: YY_k, air_mass, PP, TT, del, psat1, psat2, psf_1, psf_2 real :: air_mass_1, air_mass_2, sum1, sum2, init_water_1, init_water_2 logical :: spot_exist=.true., lmake_spot, lline_profile=.false. real :: Rgas_loc=8.314472688702992E+7, T_tmp real :: aa0= 6.107799961, aa1= 4.436518521e-1 real :: aa2= 1.428945805e-2, aa3= 2.650648471e-4 real :: aa4= 3.031240396e-6, aa5= 2.034080948e-8, aa6= 6.136820929e-11 intent(in) :: air_mass ! ! Reinitialization of T, water => rho ! if (linit_temperature) then del=(init_x2-init_x1)*0.2 if (lcurved_xz) then do j=n1,n2 init_x1_arz(j)=init_x1*(1-0.6*sin(4.*PI*z(j)/Lxyz(3))) init_x2_arz(j)=init_x2*(1+0.6*sin(4.*PI*z(j)/Lxyz(3))) enddo del_ar1z(:)=del*(1-0.6*sin(4.*PI*z(:)/Lxyz(3))) del_ar2z(:)=del*(1+0.6*sin(4.*PI*z(:)/Lxyz(3))) elseif (lcurved_xy) then do j=m1,m2 init_x1_ary(j)=init_x1*(1-0.1*sin(4.*PI*y(j)/Lxyz(2))) init_x2_ary(j)=init_x2*(1+0.1*sin(4.*PI*y(j)/Lxyz(2))) enddo del_ar1y(:)=del*(1-0.1*sin(4.*PI*y(:)/Lxyz(2))) del_ar2y(:)=del*(1+0.1*sin(4.*PI*y(:)/Lxyz(2))) else init_x1_ary=init_x1 init_x2_ary=init_x2 init_x1_arz=init_x1 init_x2_arz=init_x2 del_ar1y(:)=del del_ar2y(:)=del del_ar1z(:)=del del_ar2z(:)=del endif ! do i=l1,l2 if (x(i)<0) then del_ary=del_ar1y del_arz=del_ar1z else del_ary=del_ar2y del_arz=del_ar2z endif ! if (ltanh_prof_xy) then do j=m1,m2 f(i,j,:,ilnTT)=log((init_TT2+init_TT1)*0.5 & +((init_TT2-init_TT1)*0.5) & *(exp(x(i)/del_ary(j))-exp(-x(i)/del_ary(j))) & /(exp(x(i)/del_ary(j))+exp(-x(i)/del_ary(j)))) enddo elseif (ltanh_prof_xz) then do j=n1,n2 f(i,:,j,ilnTT)=log((init_TT2+init_TT1)*0.5 & +((init_TT2-init_TT1)*0.5) & *(1.-exp(-2.*x(i)/del_arz(j))) & /(1.+exp(-2.*x(i)/del_arz(j)))) enddo else do j=m1,m2 if (x(i)<=init_x1_ary(j)) then f(i,j,:,ilnTT)=alog(init_TT1) endif if (x(i)>=init_x2_ary(j)) then f(i,j,:,ilnTT)=alog(init_TT2) endif if (x(i)>init_x1_ary(j) .and. x(i)1) .and. (nygrid==1).and. (nzgrid==1)) then f(:,:,:,iux)=f(:,:,:,iux)+init_ux endif ! if (lroot) print*, ' Saturation Pressure, Pa ', psf_1, psf_2 if (lroot) print*, ' psf, Pa ', psf_1, psf_2 if (lroot) print*, ' pw1, Pa ', (exp(f(l1,m1,n1,ilnrho))*Rgas_loc*init_TT1/18.), PP*air_mass_1/18. if (lroot) print*, ' pw2, Pa ', (exp(f(l2,m1,n1,ilnrho))*Rgas_loc*init_TT2/18.), PP*air_mass_2/18. if (lroot) print*, ' saturated water mass fraction', psat1/PP, psat2/PP ! if (lroot) print*, 'New Air density, g/cm^3:' ! if (lroot) print '(E10.3)', PP/(k_B_cgs/m_u_cgs)*maxval(air_mass_ar)/TT if (lroot) print*, 'New Air mean weight, g/mol', maxval(air_mass_ar) if ((lroot) .and. (nx >1 )) then print*, 'density', exp(f(l1,4,4,ilnrho)), exp(f(l2,4,4,ilnrho)) endif if ((lroot) .and. (nx >1 )) then print*, 'temperature', exp(f(l1,4,4,ilnTT)), exp(f(l2,4,4,ilnTT)) endif ! endsubroutine reinitialization !********************************************************************************************************* subroutine spot_init(f,init_TT_2,init_H2O_2) ! ! Initialization of the dust spot positions and dust distribution ! ! 10-may-10/Natalia: coded ! use General, only: random_number_wrapper ! real, dimension (mx,my,mz,mfarray) :: f integer :: k, j, j1,j2,j3, lx=0,ly=0,lz=0 real :: RR, init_TT_2, init_H2O_2 real, dimension (3,spot_number) :: spot_posit logical :: spot_exist=.true. ! spot_posit(:,:)=0.0 do j=1,spot_number spot_exist=.true. lx=0;ly=0; lz=0 if (nxgrid/=1) then lx=1 call random_number_wrapper(spot_posit(1,j)) spot_posit(1,j)=spot_posit(1,j)*Lxyz(1) if ((spot_posit(1,j)-1.5*spot_sizexyz0(1)+Lxyz(1))) & spot_exist=.false. print*,'positx',spot_posit(1,j),spot_exist ! if ((spot_posit(1,j)-1.5*spot_sizexyz0(2)+Lxyz(2))) & spot_exist=.false. print*,'posity',spot_posit(2,j),spot_exist endif if (nzgrid/=1) then lz=1 call random_number_wrapper(spot_posit(3,j)) spot_posit(3,j)=spot_posit(3,j)*Lxyz(3) if ((spot_posit(3,j)-1.5*spot_sizexyz0(3)+Lxyz(3))) & spot_exist=.false. endif do j1=1,mx; do j2=1,my; do j3=1,mz RR= (lx*x(j1)-spot_posit(1,j))**2 & +ly*(y(j2)-spot_posit(2,j))**2 & +lz*(z(j3)-spot_posit(3,j))**2 RR=sqrt(RR) ! if ((RR