From a5a0480b8383d65d2d818c0e5b42c5dbf40611d0 Mon Sep 17 00:00:00 2001 From: Thomas Gastine Date: Thu, 12 Sep 2024 12:28:34 +0200 Subject: [PATCH] better handling of movie types in out_dtB_frame.f90 --- src/out_dtB_frame.f90 | 91 +++++++++++++++---------------------------- src/shtransforms.f90 | 74 +++++++++++++++++------------------ 2 files changed, 68 insertions(+), 97 deletions(-) diff --git a/src/out_dtB_frame.f90 b/src/out_dtB_frame.f90 index ab347171..c08f8288 100644 --- a/src/out_dtB_frame.f90 +++ b/src/out_dtB_frame.f90 @@ -25,10 +25,9 @@ module out_dtB_frame & PadvLM_Rloc, PadvLM_LMloc, TstrLM_Rloc, TstrLM_LMloc, & & TadvLM_Rloc, TadvLM_LMloc, TdifLM_LMloc, TomeLM_Rloc, & & TomeLM_LMloc - use movie_data, only: n_movie_type, n_movie_fields, n_movie_fields_ic, & - & n_movie_file, n_movie_const, n_movie_surface, & - & movie_const, n_movie_field_type, frames, n_movies, & - & n_movie_field_start + use movie_data, only: n_movie_fields, n_movie_fields_ic, n_movie_file, & + & n_movie_const, n_movie_surface, movie_const, & + & n_movie_field_type, frames, n_movies, n_movie_field_start use logic, only: l_cond_ic use constants, only: zero, one use radial_der_even, only: get_drNS_even @@ -50,7 +49,7 @@ subroutine calc_dtB_frame() ! !-- Local variables: - integer :: n_type, n_out, n_movie + integer :: n_out, n_movie integer :: n_fields, n_field, n_const, n_r, n_store_last integer :: n_theta, n_theta_cal, n_phi, n_o, n_or integer :: lm, l, m @@ -65,21 +64,18 @@ subroutine calc_dtB_frame() do n_movie=1,n_movies - n_type =n_movie_type(n_movie) n_fields =n_movie_fields(n_movie) n_out =n_movie_file(n_movie) n_const =n_movie_const(n_movie) n_surface=n_movie_surface(n_movie) const =movie_const(n_movie) - !--- Axisymmetric dtFL or dtAB: - if ( n_type == 31 .or. n_type == 32 .or. n_type == 33 .or. & - & n_type == 41 .or. n_type == 42 .or. n_type == 43 ) then + do n_field=1,n_fields - 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 - 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_surface == 3 ) then ! Axisymmetric movies if ( n_field_type == 20 ) then call lo2r_one%transp_lm2r(PstrLM_LMloc,PstrLM_Rloc) @@ -93,6 +89,8 @@ subroutine calc_dtB_frame() call lo2r_one%transp_lm2r(TadvLM_LMloc,TadvLM_Rloc) else if ( n_field_type == 26 ) then call lo2r_one%transp_lm2r(TdifLM_LMloc,work_Rloc) + else + cycle ! Leave the n_field loop end if do n_r=nRstart,nRstop @@ -120,21 +118,7 @@ subroutine calc_dtB_frame() end do - end do - - !--- dtBr: - else if ( n_type == 103 .or. n_type == 51 .or. n_type == 52 .or. & - & n_type == 53 .or. n_type == 91 .or. n_type == 92 .or. & - & n_type == 71 .or. n_type == 72 .or. n_type == 73 .or. & - & n_type == 80 .or. n_type == 81 .or. n_type == 82 .or. & - & n_type == 83 .or. n_type == 81 .or. n_type == 82 .or. & - & n_type == 93 .or. n_type == 60 .or. n_type == 61 .or. & - & n_type == 62 .or. n_type == 63 ) then ! non-axisymmetric fields - - 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 + else if ( n_surface == 0 .or. n_surface == 1 ) then ! 3D or rcut movies if ( n_field_type == 24 ) then ! This reduces omega effect to field production of axisymm. toroidal field: @@ -342,6 +326,8 @@ subroutine calc_dtB_frame() !--- Fieldlines for theta=const. else if ( n_field_type == 52 ) then call get_Btor(b_Rloc(:,n_r),dtBr,dtBt,rMov,.false.) + else + cycle ! Do not store anything in frames(*) end if !--- Now store the stuff to frames(*): @@ -372,11 +358,11 @@ subroutine calc_dtB_frame() end do ! r loop - end do ! LOOP over fields in movie file + end if ! which surface ? - end if + end do ! LOOP over fields in movie file - end do + end do ! Loop over movies end subroutine calc_dtB_frame !---------------------------------------------------------------------------- @@ -387,7 +373,7 @@ subroutine calc_dtB_frame_IC() ! !-- Local variables: - integer :: n_type, n_out, n_fields_oc, n_fields_ic, n_movie + integer :: n_out, n_fields_oc, n_fields_ic, n_movie integer :: n_fields, n_field, n_const, n_r, n_store_last integer :: n_theta, n_theta_cal, n_phi, n_o, n_or integer :: n_surface, n_field_type @@ -404,7 +390,6 @@ subroutine calc_dtB_frame_IC() do n_movie=1,n_movies - n_type =n_movie_type(n_movie) n_fields_oc =n_movie_fields(n_movie) n_fields_ic =n_movie_fields_ic(n_movie) n_fields =n_fields_oc+n_fields_ic @@ -413,14 +398,14 @@ subroutine calc_dtB_frame_IC() n_surface =n_movie_surface(n_movie) const =movie_const(n_movie) - !--- Axisymmetric dtFL or dtAB: - if ( n_type == 31 .or. n_type == 32 .or. n_type == 33 .or. & - & n_type == 41 .or. n_type == 42 .or. n_type == 43 ) then + do n_field=n_fields_oc+1,n_fields - do n_field=n_fields_oc+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 + + !--- Axisymmetric dtFL or dtAB: + if ( n_surface == 3 ) then - n_field_type=n_movie_field_type(n_field,n_movie) - n_store_last=n_movie_field_start(n_field,n_movie)-1 do n_r=1,n_r_ic_max @@ -433,13 +418,12 @@ subroutine calc_dtB_frame_IC() else if ( n_field_type == 26 ) then call get_dtB(dtB,TdifLMIC_LMloc,llm,ulm,1,n_r_ic_max,n_r,.true.) else - dtB(:)=0.0_cp + cycle end if if ( rank == 0 ) then n_or=n_store_last+(n_r-1)*n_theta_max - if ( n_r==1 ) print*, maxval(abs(dtB)) !--- Store in frame field: do n_theta_cal=1,n_theta_max n_theta=n_theta_cal2ord(n_theta_cal) @@ -449,21 +433,8 @@ subroutine calc_dtB_frame_IC() end do - end do - - !--- dtBr: - else if ( n_type == 103 .or. n_type == 51 .or. n_type == 52 .or. & - & n_type == 53 .or. n_type == 91 .or. n_type == 92 .or. & - & n_type == 71 .or. n_type == 72 .or. n_type == 73 .or. & - & n_type == 80 .or. n_type == 81 .or. n_type == 82 .or. & - & n_type == 83 .or. n_type == 81 .or. n_type == 82 .or. & - & n_type == 93 .or. n_type == 60 .or. n_type == 61 .or. & - & n_type == 62 .or. n_type == 63 ) then ! non-axisymmetric fields - - do n_field=n_fields_oc+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 + !--- dtBr: + else if ( n_surface == 0 .or. n_surface == 1 ) then l_loop=.true. if ( n_surface == 1 .and. const >= r_icb ) l_loop= .false. @@ -578,6 +549,8 @@ subroutine calc_dtB_frame_IC() else if ( n_field_type == 52 ) then call gather_from_lo_to_rank0(b_ic_LMloc(:,n_r), workA_1d) call get_Btor(workA_1d,dtBr,dtBt,rMov,.true.) + else + cycle end if !--- Now rank=0 store the stuff for the inner core. @@ -602,14 +575,12 @@ subroutine calc_dtB_frame_IC() end if ! l_loop ? - !--- Write frame field: - !write(n_out) (real(dtBrframe(n),kind=outp),n=1,n_field_size) - end do ! LOOP over fields in movie file + end if ! n_surface ? - end if + end do ! LOOP over fields in movie file - end do + end do ! Loop over movies end subroutine calc_dtB_frame_IC !---------------------------------------------------------------------------- diff --git a/src/shtransforms.f90 b/src/shtransforms.f90 index f76ad666..4c362513 100644 --- a/src/shtransforms.f90 +++ b/src/shtransforms.f90 @@ -102,29 +102,29 @@ end subroutine finalize_transforms !------------------------------------------------------------------------------ subroutine native_qst_to_spat(Qlm, Slm, Tlm, brc, btc, bpc, lcut) ! - ! Vector spherical harmonic transform: take Q,S,T and transform them + ! Vector spherical harmonic transform: take Q,S,T and transform them ! to a vector field ! - + !-- Input variables: complex(cp), intent(in) :: Qlm(lm_max) ! Poloidal complex(cp), intent(in) :: Slm(lm_max) ! Spheroidal complex(cp), intent(in) :: Tlm(lm_max) ! Toroidal integer, intent(in) :: lcut - + !-- Output: field on grid (theta,m) for the radial grid point nR ! and equatorially symmetric and antisymmetric contribution real(cp), intent(out) :: brc(n_theta_max,n_phi_max) real(cp), intent(out) :: btc(n_theta_max,n_phi_max) real(cp), intent(out) :: bpc(n_theta_max,n_phi_max) - - !------ Legendre Polynomials + + !------ Legendre Polynomials real(cp) :: PlmG(lm_max),PlmC(lm_max) complex(cp) :: bhG(lm_max),bhC(lm_max) complex(cp) :: tmpr(n_theta_max,n_phi_max/2+1) complex(cp) :: tmpt(n_theta_max,n_phi_max/2+1) complex(cp) :: tmpp(n_theta_max,n_phi_max/2+1) - + !-- Local variables: logical :: l_Odd complex(cp) :: brES,brEA @@ -133,7 +133,7 @@ subroutine native_qst_to_spat(Qlm, Slm, Tlm, brc, btc, bpc, lcut) real(cp) :: dm complex(cp) :: bhN1M,bhN2M,bhN,bhN1,bhN2 complex(cp) :: bhS1M,bhS2M,bhS,bhS1,bhS2 - + bhG(1)=zero bhC(1)=zero do lm=2,lm_max @@ -168,7 +168,7 @@ subroutine native_qst_to_spat(Qlm, Slm, Tlm, brc, btc, bpc, lcut) do nThetaN=nThStart,nThStop,2 ! Loop over thetas for north HS nThetaS =nThetaN+1 ! same theta but for southern HS nThetaNHS=nThetaNHS+1 ! theta-index of northern hemisph. point - + brES =zero brEA =zero !--- 6 add/mult, 26 dble words @@ -187,7 +187,7 @@ subroutine native_qst_to_spat(Qlm, Slm, Tlm, brc, btc, bpc, lcut) end if tmpr(nThetaN,mc)=brES+brEA tmpr(nThetaS,mc)=brES-brEA - + bhN1=zero bhS1=zero bhN2=zero @@ -209,7 +209,7 @@ subroutine native_qst_to_spat(Qlm, Slm, Tlm, brc, btc, bpc, lcut) bhS1M=half*bhS1 bhN2M=half*bhN2 bhS2M=half*bhS2 - + !--- 6 add/mult, 20 dble words tmpt(nThetaN,mc)=bhN1M+bhN2M bhN =bhN1M-bhN2M @@ -219,7 +219,7 @@ subroutine native_qst_to_spat(Qlm, Slm, Tlm, brc, btc, bpc, lcut) tmpp(nThetaS,mc)=-ci*bhS end do ! End global loop over nTheta end do - + !-- Zero out terms with index mc > n_m_max: if ( n_m_max < n_phi_max/2+1 ) then do mc=n_m_max+1,n_phi_max/2+1 @@ -241,7 +241,7 @@ subroutine native_qst_to_spat(Qlm, Slm, Tlm, brc, btc, bpc, lcut) btc(:,1)=real(tmpt(:,1)) bpc(:,1)=real(tmpp(:,1)) end if - + end subroutine native_qst_to_spat !------------------------------------------------------------------------------ subroutine native_sphtor_to_spat(Slm, Tlm, btc, bpc, lcut) @@ -249,23 +249,23 @@ subroutine native_sphtor_to_spat(Slm, Tlm, btc, bpc, lcut) ! Use spheroidal and toroidal potentials to transform them to angular ! vector components btheta and bphi ! - + !-- Input variables: complex(cp), intent(in) :: Slm(lm_max) complex(cp), intent(in) :: Tlm(lm_max) integer, intent(in) :: lcut - + !-- Output: field on grid (theta,m) for the radial grid point nR ! and equatorially symmetric and antisymmetric contribution real(cp), intent(out) :: btc(n_theta_max,n_phi_max) real(cp), intent(out) :: bpc(n_theta_max,n_phi_max) - - !------ Legendre Polynomials + + !------ Legendre Polynomials real(cp) :: PlmG(lm_max),PlmC(lm_max) complex(cp) :: bhG(lm_max),bhC(lm_max) complex(cp) :: tmpt(n_theta_max,n_phi_max/2+1) complex(cp) :: tmpp(n_theta_max,n_phi_max/2+1) - + !-- Local variables: logical :: l_Odd integer :: nThetaN,nThetaS,nThetaNHS @@ -273,7 +273,7 @@ subroutine native_sphtor_to_spat(Slm, Tlm, btc, bpc, lcut) real(cp) :: dm complex(cp) :: bhN1M,bhN2M,bhN,bhN1,bhN2 complex(cp) :: bhS1M,bhS2M,bhS,bhS1,bhS2 - + do lm=1,lm_max bhG(lm)=Slm(lm)-ci*Tlm(lm) bhC(lm)=Slm(lm)+ci*Tlm(lm) @@ -305,7 +305,7 @@ subroutine native_sphtor_to_spat(Slm, Tlm, btc, bpc, lcut) do nThetaN=nThStart,nThStop,2 ! Loop over thetas for north HS nThetaS =nThetaN+1 ! same theta but for southern HS nThetaNHS=nThetaNHS+1 ! theta-index of northern hemisph. point - + !--- 6 add/mult, 26 dble words do lm=lStart(mc),lmS-1,2 PlmG(lm) =dPlm(lm,nThetaNHS) -dm*Plm(lm,nThetaNHS) @@ -317,7 +317,7 @@ subroutine native_sphtor_to_spat(Slm, Tlm, btc, bpc, lcut) PlmG(lmS)=dPlm(lmS,nThetaNHS)-dm*Plm(lmS,nThetaNHS) PlmC(lmS)=dPlm(lmS,nThetaNHS)+dm*Plm(lmS,nThetaNHS) end if - + bhN1=zero bhS1=zero bhN2=zero @@ -339,7 +339,7 @@ subroutine native_sphtor_to_spat(Slm, Tlm, btc, bpc, lcut) bhS1M=half*bhS1 bhN2M=half*bhN2 bhS2M=half*bhS2 - + !--- 6 add/mult, 20 dble words tmpt(nThetaN,mc)=bhN1M+bhN2M bhN =bhN1M-bhN2M @@ -349,7 +349,7 @@ subroutine native_sphtor_to_spat(Slm, Tlm, btc, bpc, lcut) tmpp(nThetaS,mc)=-ci*bhS end do end do ! End global loop over mc - + !-- Zero out terms with index mc > n_m_max: if ( n_m_max < n_phi_max/2+1 ) then do mc=n_m_max+1,n_phi_max/2+1 @@ -368,14 +368,14 @@ subroutine native_sphtor_to_spat(Slm, Tlm, btc, bpc, lcut) btc(:,1)=real(tmpt(:,1)) bpc(:,1)=real(tmpp(:,1)) end if - + end subroutine native_sphtor_to_spat !------------------------------------------------------------------------------ subroutine native_axi_to_spat(Slm, sc) !-- Input variable complex(cp), intent(in) :: Slm(l_max+1) - + !-- Output variable real(cp), intent(out) :: sc(n_theta_max) @@ -407,26 +407,26 @@ subroutine native_toraxi_to_spat(Tlm, btc, bpc) ! Use spheroidal and toroidal potentials to transform them to angular ! vector components btheta and bphi ! - + !-- Input variables: complex(cp), intent(in) :: Tlm(l_max+1) - + !-- Output: field on grid (theta,m) for the radial grid point nR ! and equatorially symmetric and antisymmetric contribution real(cp), intent(out) :: btc(n_theta_max) real(cp), intent(out) :: bpc(n_theta_max) - - !------ Legendre Polynomials + + !------ Legendre Polynomials real(cp) :: PlmG(l_max+1),PlmC(l_max+1) complex(cp) :: bhG(l_max+1),bhC(l_max+1) - + !-- Local variables: integer :: nThetaN,nThetaS,nThetaNHS integer :: lm,lmS,l - + complex(cp) :: bhN1M,bhN2M,bhN,bhN1,bhN2 complex(cp) :: bhS1M,bhS2M,bhS,bhS1,bhS2 - + do l=1,l_max+1 bhG(l)=-ci*Tlm(l) bhC(l)= ci*Tlm(l) @@ -436,7 +436,7 @@ subroutine native_toraxi_to_spat(Tlm, btc, bpc) do nThetaN=1,n_theta_max,2 ! Loop over thetas for north HS nThetaS =nThetaN+1 ! same theta but for southern HS nThetaNHS=nThetaNHS+1 ! theta-index of northern hemisph. point - + lmS=lStop(1) !--- 6 add/mult, 26 dble words do lm=lStart(1),lmS-1,2 @@ -449,7 +449,7 @@ subroutine native_toraxi_to_spat(Tlm, btc, bpc) PlmG(lmS)=dPlm(lmS,nThetaNHS) PlmC(lmS)=dPlm(lmS,nThetaNHS) end if - + lmS=lStop(1) bhN1=zero bhS1=zero @@ -481,7 +481,7 @@ subroutine native_toraxi_to_spat(Tlm, btc, bpc) bpc(nThetaN)=real(-ci*bhN) bpc(nThetaS)=real(-ci*bhS) end do - + end subroutine native_toraxi_to_spat !------------------------------------------------------------------------------ subroutine native_sph_to_spat(Slm, sc, lcut) @@ -542,7 +542,7 @@ subroutine native_sph_to_spat(Slm, sc, lcut) do nThetaN=nThstart,nThStop tmp(nThetaN,mc)=zero end do ! loop over nThetaN (theta) - end do + end do end if !$omp end parallel @@ -717,7 +717,7 @@ subroutine native_spat_to_sph(scal,f1LM,lcut) f1LM(lmS)=f1LM(lmS) + f1ES1*wPlm(lmS,nTheta1) + f1ES2*wPlm(lmS,nTheta2) end if end do - end do + end do !$omp end parallel end subroutine native_spat_to_sph @@ -853,7 +853,7 @@ subroutine native_spat_to_sph_tor(vt,vp,f1LM,f2LM,lcut) end if end do ! loop over theta in block - end do + end do !$omp end parallel !-- Division by l(l+1) except for (l=0,m=0)