diff --git a/src/movie.f90 b/src/movie.f90 index dbbc6747..72c18a61 100644 --- a/src/movie.f90 +++ b/src/movie.f90 @@ -10,7 +10,7 @@ module movie_data use radial_data, only: nRstart,nRstop, n_r_icb, n_r_cmb, radial_balance use radial_functions, only: r_cmb, r_icb, r, r_ic use horizontal_data, only: theta_ord, phi, n_theta_ord2cal - use output_data, only: n_log_file, log_file, tag + use output_data, only: n_log_file, log_file, tag, n_s_max use charmanip, only: capitalize, delete_string, dble2str use useful, only: logWrite, abortRun use constants, only: pi, one @@ -30,8 +30,7 @@ module movie_data character(len=80), public :: movie(n_movies_max) ! Only for input character(len=72), public :: movie_file(n_movies_max) - logical, public :: lStoreMov(n_movies_max),lICField(n_movies_max) - logical, public :: lGeosField(n_movies_max) + logical, public :: lICField(n_movies_max), lGeosField(n_movies_max) logical :: lAxiField(n_movies_max), lPhaseField(n_movies_max) integer, public :: n_movies integer, public :: n_movie_surface(n_movies_max) @@ -336,7 +335,7 @@ subroutine get_movie_type integer :: n_fields,n_fields_oc,n_fields_ic integer :: n_field_size, n_field_size_ic, n_field_start integer :: n_field_type(n_movie_fields_max) - logical :: lStore, lIC, lAxi, lGeos, lPhase + logical :: lIC, lAxi, lGeos, lPhase !--- Initialize first storage index: n_field_start=1 @@ -358,7 +357,6 @@ subroutine get_movie_type do i=1,n_movies_max - lStore=.true. lIC =.false. lPhase=.false. lGeos =.false. @@ -688,7 +686,6 @@ subroutine get_movie_type typeStr=' geos phi-component of velocity ' file_name='geosVPHI_' n_fields=1 - lStore=.false. l_geosMovie=.true. n_field_type(1)=101 lGeos=.true. @@ -734,7 +731,6 @@ subroutine get_movie_type typeStr=' geos z-component of vorticity ' file_name='geosVorZ_' n_fields=1 - lStore=.false. l_geosMovie=.true. n_field_type(1)=102 lGeos=.true. @@ -774,7 +770,6 @@ subroutine get_movie_type typeStr=' geos s-component of velocity ' file_name='geosVS_' n_fields=1 - lStore=.false. l_geosMovie=.true. n_field_type(1)=100 lGeos=.true. @@ -969,7 +964,7 @@ subroutine get_movie_type else if ( lGeos ) then ! Geos average n_surface=-2 ! constant theta n_const=1 ! - n_field_size=n_phi_max*n_r_max + n_field_size=n_phi_max*n_s_max n_field_size_ic=0 const=r_cmb else if ( lPhase ) then ! Melting radius or temp gradient at rm @@ -1132,7 +1127,6 @@ subroutine get_movie_type !--- Now store the necessary information: !------ Increase number of movies: n_movies=n_movies+1 - lStoreMov(n_movies)=lStore lICField(n_movies)=lIC lGeosField(n_movies)=lGeos lAxiField(n_movies)=lAxi @@ -1208,27 +1202,17 @@ subroutine get_movie_type do n=1,n_fields if ( n_fields_oc > 0 ) then n_movie_field_type(n,n_movies) =n_field_type(n) - if ( lStore ) then - n_movie_field_start(n,n_movies)=n_field_start - n_field_start=n_field_start+n_field_size - n_movie_field_stop(n,n_movies) =n_field_start-1 - l_store_frame=.true. - else - n_movie_field_start(n,n_movies)=-1 - n_movie_field_stop(n,n_movies) =-1 - end if + n_movie_field_start(n,n_movies)=n_field_start + n_field_start=n_field_start+n_field_size + n_movie_field_stop(n,n_movies) =n_field_start-1 + l_store_frame=.true. end if if ( n_fields_ic > 0 ) then n_ic=n_fields_oc+n n_movie_field_type(n_ic,n_movies)=n_field_type(n) - if ( lStore ) then - n_movie_field_start(n_ic,n_movies)=n_field_start - n_field_start=n_field_start+n_field_size_ic - n_movie_field_stop(n_ic,n_movies)=n_field_start-1 - else - n_movie_field_start(n_ic,n_movies)=-1 - n_movie_field_stop(n_ic,n_movies) =-1 - end if + n_movie_field_start(n_ic,n_movies)=n_field_start + n_field_start=n_field_start+n_field_size_ic + n_movie_field_stop(n_ic,n_movies)=n_field_start-1 end if end do @@ -1308,6 +1292,10 @@ subroutine movie_gather_frames_to_rank0() select case(n_surface) + case(-2) ! Cylindrical geos movie + ! Data is already on rank 0 + cycle + case(-1) ! Earth Surface ! theta-phi surface for n_r=1 (CMB) ! frames is already existent on rank 0 with all diff --git a/src/outGeos.f90 b/src/outGeos.f90 index e636cdd6..8df7e334 100644 --- a/src/outGeos.f90 +++ b/src/outGeos.f90 @@ -14,10 +14,11 @@ module geos use mem_alloc, only: bytes_allocated use radial_data, only: radial_balance, nRstart, nRstop use radial_functions, only: or1, or2, r_ICB, r_CMB, r, orho1, orho2, beta - use output_data, only: sDens, zDens, tag + use output_data, only: sDens, zDens, tag, n_s_max use horizontal_data, only: n_theta_cal2ord, O_sin_theta_E2, theta_ord, & & O_sin_theta, cosTheta, sinTheta - use movie_data, only: n_movie_field_type, n_movie_fields, n_movie_file + use movie_data, only: n_movie_field_type, n_movie_fields, n_movie_file, & + & n_movies, n_movie_field_start, frames use truncation, only: n_phi_max, n_theta_max, n_r_max, nlat_padded, l_max use integration, only: simps, cylmean_otc, cylmean_itc use logic, only: l_save_out @@ -38,12 +39,11 @@ module geos integer :: nPstop ! Stoping nPhi index when MPI distributed type(load), allocatable :: phi_balance(:) ! phi-distributed balance integer :: n_geos_file ! file unit for geos.TAG - integer, public :: n_s_max ! Number of cylindrical points integer :: n_s_otc ! Index for last point outside TC character(len=72) :: geos_file ! file name public :: initialize_geos, finalize_geos, calcGeos, outGeos, outOmega, & - & write_geos_frame + & calc_geos_frame contains @@ -84,8 +84,6 @@ subroutine initialize_geos(l_geos, l_SRIC, l_geosMovie) if ( l_geos .or. l_SRIC .or. l_geosMovie ) then !-- Cylindrical radius - n_s_max = n_r_max+int(r_ICB*n_r_max) - n_s_max = int(sDens*n_s_max) allocate( cyl(n_s_max) ) bytes_allocated=bytes_allocated+n_s_max*SIZEOF_DEF_REAL @@ -494,57 +492,67 @@ subroutine outGeos(time,Geos,GeosA,GeosZ,GeosM,GeosNAP,Ekin) end subroutine outGeos !------------------------------------------------------------------------------------ - subroutine write_geos_frame(n_movie) + subroutine calc_geos_frame() ! ! This subroutine handles the computation and the writing of geos movie ! files. ! - !-- Input variables - integer, intent(in) :: n_movie ! The index of the movie in list of movies - !-- Local variables - integer :: n_p, n_fields, n_field, n_field_type, n_out + integer :: n_p, n_fields, n_field, n_field_type, n_out, n_movie + integer :: n_store_last, n_phi, n_s, n_o real(cp) :: dat(n_s_max,nPstart:nPstop), dat_full(n_phi_max,n_s_max) real(cp) :: datITC_N(n_s_max), datITC_S(n_s_max) - n_fields=n_movie_fields(n_movie) - n_out =n_movie_file(n_movie) + do n_movie=1,n_movies + n_fields=n_movie_fields(n_movie) + n_out =n_movie_file(n_movie) - do n_field=1,n_fields - n_field_type=n_movie_field_type(n_field,n_movie) + do n_field=1,n_fields + n_field_type=n_movie_field_type(n_field,n_movie) + n_store_last=n_movie_field_start(n_field,n_movie)-1 - if ( n_field_type == 100 ) then ! Us - call transp_R2Phi(us_Rloc, us_Ploc) - !-- Z averaging - do n_p=nPstart,nPstop - call cylmean(us_Ploc(:,:,n_p),dat(:,n_p),datITC_N,datITC_S) - dat(:,n_p)=dat(:,n_p)+datITC_N(:) - end do - else if (n_field_type == 101 ) then ! Uphi - call transp_R2Phi(up_Rloc, up_Ploc) - !-- Z averaging - do n_p=nPstart,nPstop - call cylmean(up_Ploc(:,:,n_p),dat(:,n_p),datITC_N,datITC_S) - dat(:,n_p)=dat(:,n_p)+datITC_N(:) - end do - else if (n_field_type == 102 ) then ! Vort z - call transp_R2Phi(wz_Rloc, wz_Ploc) - !-- Z averaging - do n_p=nPstart,nPstop - call cylmean(wz_Ploc(:,:,n_p),dat(:,n_p),datITC_N,datITC_S) - dat(:,n_p)=dat(:,n_p)+datITC_N(:) - end do - end if + if ( n_field_type == 100 ) then ! Us + call transp_R2Phi(us_Rloc, us_Ploc) + !-- Z averaging + do n_p=nPstart,nPstop + call cylmean(us_Ploc(:,:,n_p),dat(:,n_p),datITC_N,datITC_S) + dat(:,n_p)=dat(:,n_p)+datITC_N(:) + end do + else if (n_field_type == 101 ) then ! Uphi + call transp_R2Phi(up_Rloc, up_Ploc) + !-- Z averaging + do n_p=nPstart,nPstop + call cylmean(up_Ploc(:,:,n_p),dat(:,n_p),datITC_N,datITC_S) + dat(:,n_p)=dat(:,n_p)+datITC_N(:) + end do + else if (n_field_type == 102 ) then ! Vort z + call transp_R2Phi(wz_Rloc, wz_Ploc) + !-- Z averaging + do n_p=nPstart,nPstop + call cylmean(wz_Ploc(:,:,n_p),dat(:,n_p),datITC_N,datITC_S) + dat(:,n_p)=dat(:,n_p)+datITC_N(:) + end do + else + cycle + end if - !-- MPI gather - call gather_Ploc(dat, dat_full) + !-- MPI gather + call gather_Ploc(dat, dat_full) - !-- Write outputs - if ( rank == 0 ) write(n_out) real(dat_full,kind=outp) - end do + !-- Write outputs + if ( rank == 0 ) then + do n_s=1,n_s_max + n_o=n_store_last+(n_s-1)*n_phi_max + do n_phi=1,n_phi_max + frames(n_phi+n_o)=real(dat_full(n_phi,n_s),kind=outp) + end do + end do + end if + end do ! loop over n_fields + end do ! loop over n_movies - end subroutine write_geos_frame + end subroutine calc_geos_frame !------------------------------------------------------------------------------------ subroutine outOmega(z, omega_ic) ! diff --git a/src/out_TO.f90 b/src/out_TO.f90 index dcf9afe2..ee9e3f41 100644 --- a/src/out_TO.f90 +++ b/src/out_TO.f90 @@ -12,7 +12,7 @@ module outTO_mod use truncation, only: n_r_max, n_theta_max, n_phi_max, minc use mem_alloc, only: bytes_allocated use num_param, only: tScale - use output_data, only: sDens, zDens, tag, runid, log_file, n_log_file + use output_data, only: sDens, zDens, tag, runid, log_file, n_log_file, n_s_max use radial_data, only: radial_balance, nRstart, nRstop use radial_functions, only: r_ICB, r_CMB, r use horizontal_data, only: theta_ord, phi @@ -39,7 +39,7 @@ module outTO_mod real(cp), allocatable :: dzLFAS(:,:), V2AS(:,:), Bs2AS(:,:), BspAS(:,:), VAS(:,:) real(cp), allocatable :: BspdAS(:,:), BpsdAS(:,:), BzpdAS(:,:), BpzdAS(:,:) real(cp), allocatable :: dzPenAS(:,:) - integer :: n_s_otc, n_s_max, n_NHS_file, n_SHS_file, n_TOmov_file, n_TO_file + integer :: n_s_otc, n_NHS_file, n_SHS_file, n_TOmov_file, n_TO_file real(cp) :: volcyl_oc character(len=64) :: movFile, TOFile integer :: nTOmovSets ! Number of TO_mov frames @@ -89,8 +89,6 @@ subroutine initialize_outTO_mod() end if !-- Cylindrical radius - n_s_max = n_r_max+int(r_ICB*n_r_max) - n_s_max = int(sDens*n_s_max) allocate( cyl(n_s_max) ) bytes_allocated=bytes_allocated+n_s_max*SIZEOF_DEF_REAL diff --git a/src/out_movie_file.f90 b/src/out_movie_file.f90 index cf315e36..f2904d3e 100644 --- a/src/out_movie_file.f90 +++ b/src/out_movie_file.f90 @@ -1,8 +1,12 @@ module out_movie + ! + ! This module handles the storage of the relevant diagnostics in the public + ! array frames(*) and then the writing of the corresponding movie files + ! use precision_mod use parallel_mod, only: rank - use geos, only: cyl, n_s_max, write_geos_frame + use geos, only: cyl use truncation, only: n_phi_max, n_theta_max, minc, lm_max, l_max, & & n_m_max, lm_maxMag, n_r_maxMag, n_r_ic_maxMag, & & n_r_ic_max, n_r_max, nlat_padded @@ -10,9 +14,8 @@ module out_movie & n_movie_const, n_movie_field_type, lGeosField, & & n_movie_field_start,n_movie_field_stop, & & movieDipColat, movieDipLon, movieDipStrength, & - & movieDipStrengthGeo, movie_const, & - & lStoreMov, n_movie_file, n_movie_fields_ic, & - & movie_file + & movieDipStrengthGeo, movie_const, movie_file, & + & n_movie_file, n_movie_fields_ic use radial_data, only: n_r_icb, n_r_cmb use radial_functions, only: orho1, orho2, or1, or2, or3, or4, beta, & & r_surface, r_cmb, r, r_ic, temp0 @@ -26,7 +29,7 @@ module out_movie use sht, only: torpol_to_spat, toraxi_to_spat use logic, only: l_save_out, l_cond_ic, l_mag, l_full_sphere use constants, only: zero, half, one, two - use output_data, only: runid + use output_data, only: runid, n_s_max use useful, only: abortRun implicit none @@ -251,23 +254,17 @@ subroutine write_movie_frame(n_frame,time,omega_ic,omega_ma) end if !------ Write frames: - if ( .not. lStoreMov(n_movie) ) then - if ( lGeosField(n_movie) ) then - call write_geos_frame(n_movie) - end if - else - if ( rank == 0 ) then - do n_field=1,n_fields - n_start=n_movie_field_start(n_field,n_movie) - if ( n_fields_oc > 0 ) then - n_stop =n_movie_field_stop(n_field,n_movie) - end if - if ( n_fields_ic > 0 ) then - n_stop=n_movie_field_stop(n_fields_oc + n_field,n_movie) - end if - write(n_out) (real(frames(n),kind=outp),n=n_start,n_stop) - end do - end if + if ( rank == 0 ) then + do n_field=1,n_fields + n_start=n_movie_field_start(n_field,n_movie) + if ( n_fields_oc > 0 ) then + n_stop =n_movie_field_stop(n_field,n_movie) + end if + if ( n_fields_ic > 0 ) then + n_stop=n_movie_field_stop(n_fields_oc + n_field,n_movie) + end if + write(n_out) (real(frames(n),kind=outp),n=n_start,n_stop) + end do end if if ( l_save_out .and. rank==0 ) close(n_out) diff --git a/src/output.f90 b/src/output.f90 index 8b2afa97..f7b3c111 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -24,7 +24,7 @@ module output_mod & l_perpPar, l_energy_modes, l_heat, l_hel, l_par, & & l_chemical_conv, l_movie, l_full_sphere, l_spec_avg, & & l_phase_field, l_hemi, l_dtBmovie, l_phaseMovie, & - & l_dtPhaseMovie + & l_dtPhaseMovie, l_geosMovie use fields, only: omega_ic, omega_ma, b_ic,db_ic, aj_ic, & & w_LMloc, dw_LMloc, ddw_LMloc, p_LMloc, xi_LMloc, & & s_LMloc, ds_LMloc, z_LMloc, dz_LMloc, b_LMloc, & @@ -44,7 +44,7 @@ module output_mod use constants, only: vol_oc, vol_ic, mass, surf_cmb, two, three, zero use outMisc_mod, only: outHeat, outHelicity, outHemi, outPhase, get_onset, & & calc_melt_frame - use geos, only: outGeos, outOmega + use geos, only: outGeos, outOmega, calc_geos_frame use outRot, only: write_rot use integration, only: rInt_R use outPar_mod, only: outPar, outPerpPar @@ -659,6 +659,8 @@ subroutine output(time,tscheme,n_time_step,l_stop_time,l_pot,l_log, & if ( l_phaseMovie .or. l_dtphaseMovie ) call calc_melt_frame() + if ( l_geosMovie ) call calc_geos_frame() + call movie_gather_frames_to_rank0() if ( l_movie_ic .and. l_store_frame ) then diff --git a/src/output_data.f90 b/src/output_data.f90 index 9acf4117..bac0404f 100644 --- a/src/output_data.f90 +++ b/src/output_data.f90 @@ -68,6 +68,7 @@ module output_data character(len=72), public :: log_file character(len=72), public :: lp_file !----- Z-integrated output: + integer, public :: n_s_max ! Number of radial grid points in s-direction real(cp), public :: sDens ! Density in s when using z-integration real(cp), public :: zDens ! Density in z when using z-integration diff --git a/src/preCalculations.f90 b/src/preCalculations.f90 index c299d084..38468bd3 100644 --- a/src/preCalculations.f90 +++ b/src/preCalculations.f90 @@ -765,6 +765,10 @@ subroutine preCalc(tscheme) call abortRun('LCR not compatible with imposed field!') end if + !-- Define n_s_max for cylindral grid + n_s_max = n_r_max+int(r_icb*n_r_max) + n_s_max = int(sDens*n_s_max) + if (l_radial_flow_bc) then if ( ellipticity_cmb /= 0.0_cp ) then ellip_fac_cmb=-two*r_cmb**3*ellipticity_cmb*omega_ma1*omegaOsz_ma1