better handling of geos movies: now make use of frames(*)

This allows to simplify the handling of movie outputs in out_movie_file.f90

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

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)
       !

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
 

src/out_movie_file.f90

@@ -1,18 +1,21 @@
 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
    use movie_data, only: frames, n_movie_fields, n_movies, n_movie_surface, &
        &                 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)

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

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