module coef_module use params_module, only: nmlon,nmlon_h,nmlat_h,nhgt_fix,nlonlat,nmlatS2_h,nmlat_T1, & bijSum,jlatm_JT, & rho,rho_s,pccolatrad,pc_uncoupled,rtd,ylatm,& ! art_2023.06.20 test_loc_pot ! local test potential is used implicit none ! real :: coef(nmlon,nmlat_h,nhgt_fix,9,2) ! for each P-point each hemisphere and height ! real :: coef_ns(nmlon,nmlat_T1,10) ! lhs+rhs: for each P-point ! real :: coef_ns2(nmlon,nmlat_h,2,10) ! coef array summed up in altitude & rhs added ! real(8) :: rhs_ns(nlonlat) ! forcing: for each P-point ! real(8) :: lhs_ns(nlonlat,nlonlat) ! lhs: for each P-point real, allocatable :: coef(:,:,:,:,:) ! for each P-point each hemisphere and height real, allocatable :: coef_ns(:,:,:) ! lhs+rhs: for each P-point real, allocatable :: coef_ns2(:,:,:,:) ! lhs+rhs: for each P-point real(8), allocatable :: rhs_ns(:) ! forcing: for each P-point real(8), allocatable :: lhs_ns(:,:) ! lhs: for each P-point real(8), allocatable :: lhs_ns_wob(:,:) ! lhs: for each P-point lhs wiithout bij values needed for test test_loc_pot =.true. real(8), allocatable :: fac_hlg(:,:) ! fac: for each P-point - this was before on fieldline_p points ! add the g for global and leave the old one in there unchanged (can be removed once it is working) real,parameter :: phi_pol = 0. ! real,parameter :: phi_pol = 1. ! Art_2023.06.20 For test5 with antisymmetric polar potentials contains !----------------------------------------------------------------------------- subroutine calc_coef ! use fieldline_s_module,only: fieldline_s1,fline_s1, & fieldline_s2,fline_s2 use fieldline_p_module,only: fieldline_p,fline_p implicit none ! ! coef ordering from tiegcm ! ^ equatorward ! coef(4) (i-1,j+1) coef(3) (i,j+1) coef(2) (i+1,j+1) ! coef(5) (i-1,j) coef(9) (i,j) coef(1) (i+1,j) ! coef(6) (i-1,j-1) coef(7) (i,j-1) coef(8) (i+1,j-1) ! v poleward ! ! relationship between P,S1, and S2 point for the same index (i,j) ! P(i,j) then is really S1(i+0.5,j) and S2(i,j+0.5) with j increasing equatorward ! coefficient is calculated at P points ! integer :: isn,i,j,k,im,nmax,status,ic,ip,ii real :: N2p_p,N2h_p ! allocate arrays allocate(coef(nmlon,nmlat_h,nhgt_fix,9,2),STAT=status) ! for each P-point each hemisphere and height if(status /= 0 ) write(6,*) 'alloc coef failed' allocate(coef_ns(nmlon,nmlat_T1,10),STAT=status) ! lhs+rhs: for each P-point 2023.02 am solve two hemisphere if(status /= 0 ) write(6,*) 'alloc coef_ns failed' allocate(rhs_ns(nlonlat),STAT=status) ! forcing: for each P-point if(status /= 0 ) write(6,*) 'alloc rhs_ns failed' allocate(lhs_ns(nlonlat,nlonlat),STAT=status) ! lhs: for each P-point if(status /= 0 ) write(6,*) 'alloc lhs_ns failed' allocate(coef_ns2(nmlon,nmlat_h,2,10),STAT=status) ! lhs+rhs: for each P-point both hemispheres if(status /= 0 ) write(6,*) 'alloc coef_ns failed' allocate(fac_hlg(nlonlat,nlonlat),STAT=status) ! fac_hlg: for each P-point if(status /= 0 ) write(6,*) 'alloc fac_hlg failed' coef = 0. ! am 01.2023_pole initilize to zero coef_ns2 = 0. do isn = 1,2 ! loop over both hemisphere do i=1,nmlon ! loop over all longitudes if(i == 1) then im = nmlon else im = i-1 endif if(i == nmlon) then ! am 01.2023_pole need also i+1 values ip = 1 else ip = i+1 endif ! am 01.2023_pole j=1 nmax = fline_p(i,j,isn)%npts ! maximum of points on fieldline ! S-pole Eq 338 Art Jan. 17, 2023 notes ! Sum_i C3(i,1) * Phi(i,j+1) + C9P * Phi^P + beta * Phi^P* + I3^TP = 0 ! note for I3^TP the sign will be changed since in the code it is on the RHS ! ! C3(i,1)* Phi(i,j+1) -> Phi(i,j+1) * Sum_k=1^K [N2P(i,3/2,k)-N2H(i+1,3/2,k)+N2H(i-1,3/2,k)] -> coef(3) S-pole ! -Sum_k=1^K[ Sum_i=1^nmlon N2P(i,3/2,k)]] - Sum_i^nmlon [b(i,1] -> coef(9) S-pole ! -Sum_k=1^K[ Sum_i=1^nmlon S(i,1,k)]] + Sum_i^nmlon [b(i,1]Phi^NP -> coef(10) S-pole ! beta = Sum_i b(i,1) -> precaluclated ! I3^TP = Sum_i [ Sum_k^K S(i,1,k) - I3^R(i,1) ] + I3^P1/2 ! with I3^P(1/2) = Sum_i I3i,1,1/2 ! I3^R(i,1) = M3(i,1,K+1/2) * Jr^R(i,1) -> upper boundary is calculated later in calc_fac ! S(i,1,k) =-M2(i,3/2,k)*Je2^D(i,3/2,k) -> calculated in calc_S and already has minus sign for RHS included ! includes I3^P(1/2) at S(i,1,1) ! ! Phi^SP(i=1,j=1) = Phi(i,j=1) for i=2,nmlon -> C9N(i,1) = 1. C*(1,1) = -1 C1S-C8S = 0. C10S=0. ! in NH this is overwritten later Phi(i,1) = Phi^NP -> C9N(i,1) = 1. C10N(i,1) = Phi^NP ! C1N-C8N = 0. ! coef(nmlon,nmlat_h,nhgt_fix,9,2) ! coef_ns2(nmlon,nmlat_h,2,10) ! if(isn.eq.1) then ! only for S-pole art_2023.06.20 art need coef_ns2 for n-pole do k=1,nmax coef_ns2(i,j,isn,3) = coef_ns2(i,j,isn,3) + fline_s2(im,j,isn)%N2h(k) + & ! need to move the C3(i,j) to the approriate place in LHS fline_s2(i,j,isn)%N2p(k)- fline_s2(ip,j,isn)%N2h(k) ! N2P(i,3/2,k)-N2H(i+1,3/2,k)+N2H(i-1,3/2,k) coef_ns2(1,j,isn,9) = coef_ns2(1,j,isn,9) - fline_s2(i,j,isn)%N2p(k) ! -Sum_k=1^K Sum_i=1^nmlon N2P(i,3/2,k)]] -> put into i=1 coef_ns2(1,j,9,isn) coef_ns2(1,j,isn,10) = coef_ns2(1,j,isn,10)+ fline_p(i,j,isn)%S(k) ! Sum_k=1^K Sum_i=1^nmlon S(i,1,k) -> put into i=1 coef_ns2(1,j,1-,isn) end do ! end height loop ! art_2023.06.20 comment this out - structural change ! Since calculation of bij needs coef_ns2, move the calculation of bij effects to subroutine const_rhs ! coef_ns2(1,j,isn,9) = coef_ns2(1,j,isn,9) - fline_p(i,j,isn)%bij ! -Sum_i^nmlon [b(i,1] ! coef_ns2(1,j,isn,10)= coef_ns2(1,j,isn,10)- fline_p(i,j,isn)%bij*phi_pol ! Sum_i^nmlon [b(i,1] * phi_pol -> minus sign since on RHS while in notes on LHS ! ! ! don't change coef_ns2 here. Values for ic=3 are needed in const_rhs. ! if(i.ne.1) then ! for i /= 1 ! coef_ns2(i,j,isn,9) = 1. ! Phi^SP(i=1,j=1) = Phi(i,j=1) for i=2,nmlon ! coef_ns2(i,j,isn,1:8) = 0. ! need to set in matrix value at entry for phi(i=1,j=1) SH ! coef_ns2(i,j,isn,10) = 0. ! endif ! else ! N-pole Also need coef_ns2 for N pole ! coef_ns2(i,j,isn,10) = phi_pol ! N-pole for each i Phi^N(i,1)=Phi^NP ! coef_ns2(i,j,isn,9) = 1 ! coef_ns2(i,j,isn,1:8)= 0 ! endif ! end of art_2023.06.20 change ! am 01.2023_pole end do j=2,nmlat_h ! loop over all latitudes in one hemisphere nOT THE POLE nmax = fline_p(i,j,isn)%npts ! maximum of points on fieldline do k=1,nmax ! longitudinal wrap around if((nmlat_h-j+1) == k) then ! top volume at equator N2h_p = 0. N2p_p = 0. else N2h_p = fline_s2(i,j,isn)%N2h(k) ! i,j+0.5 N2p_p = fline_s2(i,j,isn)%N2p(k) ! i,j+0.5 endif coef(i,j,k,2,isn) = -fline_s1(i,j,isn)%N1h(k) + N2h_p coef(i,j,k,3,isn) = fline_s1(im,j,isn)%N1h(k) - & fline_s1(i,j,isn)%N1h(k)+ N2p_p coef(i,j,k,4,isn) = fline_s1(im,j,isn)%N1h(k) - N2h_p ! coef(i,j,k,5,isn) = fline_s1(im,j,isn)%N1p(k) + & fline_s2(i,j-1,isn)%N2h(k)-N2h_p coef(i,j,k,6,isn) = -fline_s1(im,j,isn)%N1h(k) + & fline_s2(i,j-1,isn)%N2h(k) coef(i,j,k,7,isn) = -fline_s1(im,j,isn)%N1h(k) + & fline_s1(i,j,isn)%N1h(k)+fline_s2(i,j-1,isn)%N2p(k) coef(i,j,k,8,isn) = fline_s1(i,j,isn)%N1h(k) - & fline_s2(i,j-1,isn)%N2h(k) coef(i,j,k,1,isn) = fline_s1(i,j,isn)%N1p(k) - & fline_s2(i,j-1,isn)%N2h(k)+N2h_p coef(i,j,k,9,isn) = -fline_s1(im,j,isn)%N1p(k) - & fline_s1(i,j,isn)%N1p(k)-fline_s2(i,j-1,isn)%N2p(k) -& N2p_p end do ! end height loop end do ! end lat/fieldline loop end do ! end longitude loop end do ! end hemisphere loop ! ! add the coefficients in height to get coefficients for each hemisphere ! needed for calculating high latitude FAC ! do i=1,nmlon ! loop over all longitudes do isn = 1,2 ! hemisphere loop do j=2,nmlat_h ! loop over all latitudes in one hemisphere NOT THE POLE done above nmax = fline_p(i,j,1)%npts ! maximum of points on fieldline do k=1,nmax ! height loop do ic = 1,9 ! 9-point stencil coef_ns2(i,j,isn,ic) = coef_ns2(i,j,isn,ic)+ coef(i,j,k,ic,isn) end do coef_ns2(i,j,isn,10) = coef_ns2(i,j,isn,10)+ fline_p(i,j,isn)%S(k) end do ! end height loop ! end do ! end lat/fieldline loop ! !Add polar values 230619 begin art_2023/06/20 moved polar value here j=1 ! sumS(isn) = sumS(isn) + coef_ns2(i,j,isn,10) ! sumabsS(isn) = sumabsS(isn) + abs(coef_ns2(i,j,isn,10)) ! Overwrite N-pole values if(isn.eq.2) then ! only for N-pole coef_ns2(i,j,isn,10) = phi_pol ! N-pole for each i Phi^N(i,1)=Phi^NP coef_ns2(i,j,isn,9) = 1 coef_ns2(i,j,isn,1:8)= 0 endif !230619 end ! do j=1,nmlat_h ! if(i.eq.13.and.isn.eq.1) write(26,'("coefS ",i4,9(x,e15.8))') j,(coef_ns2(i,j,isn,ii),ii=1,9) ! if(i.eq.13.and.isn.eq.2) write(27,'("coefN ",i4,9(x,e15.8))') j,(coef_ns2(i,j,isn,ii),ii=1,9) ! end do ! end do ! end hemisphere loop ! end do ! end longitude loop ! end subroutine calc_coef !----------------------------------------------------------------------------- subroutine set_coef_ns ! set the coefficients in the coefficient matrix in both hemispheres ! decide where to add the SH and NH stnecil and where not. ! change direction of stencil in the NH from coef_ns2 to coef_ns ! output: coef_ns(nmlon,nmlat_T1,10) ! lhs+rhs: for each P-point 2023.02 am solve two hemisphere use fieldline_p_module,only: fieldline_p,fline_p implicit none ! integer :: i,j,k,ic,nmax,status,jTS,jTN ! ! set the coefficient matrix and set equatorial boundary condition do i=1,nmlon ! loop over all longitudes do j=1,jlatm_JT-1 ! loop from pole to latm_JT jTS = j ! overall index SH is the same since from S-pole to equator jTN = nmlat_T1-j+1 ! overall index NH from N-pole to equator do ic = 1,9 ! 9-point stencil isn=1 SH; isn=2 NH coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic) ! SH coef_ns(i,jTN,ic) = coef_ns2(i,j,2,ic) ! NH end do ! coef_ns(i,jTN,1) = coef_ns2(i,j,2,1) ! NH 02.23.2023 ! coef_ns(i,jTN,8) = coef_ns2(i,j,2,8) ! 2 ! coef_ns(i,jTN,7) = coef_ns2(i,j,2,7) ! 3 ! coef_ns(i,jTN,6) = coef_ns2(i,j,2,6) ! 4 ! coef_ns(i,jTN,5) = coef_ns2(i,j,2,5) ! coef_ns(i,jTN,4) = coef_ns2(i,j,2,4) ! 6 ! coef_ns(i,jTN,3) = coef_ns2(i,j,2,3) ! 7 ! coef_ns(i,jTN,2) = coef_ns2(i,j,2,2) ! 8 ! coef_ns(i,jTN,9) = coef_ns2(i,j,2,9) ! coef_ns(i,jTS,10)= coef_ns2(i,j,1,10) ! SH coef_ns(i,jTN,10)= coef_ns2(i,j,2,10) ! NH ! end do ! end lat/fieldline loop j = jlatm_JT jTS = j ! overall index SH is the same since from S-pole to equator jTN = nmlat_T1-j+1 ! overall index NH from N-pole to equator coef_ns(i,jTS,5) = coef_ns2(i,j,1,5) + coef_ns2(i,j,2,5) ! added values SH + NH coef_ns(i,jTS,4) = coef_ns2(i,j,1,4) + coef_ns2(i,j,2,4) coef_ns(i,jTS,3) = coef_ns2(i,j,1,3) + coef_ns2(i,j,2,3) coef_ns(i,jTS,2) = coef_ns2(i,j,1,2) + coef_ns2(i,j,2,2) coef_ns(i,jTS,1) = coef_ns2(i,j,1,1) + coef_ns2(i,j,2,1) coef_ns(i,jTS,9) = coef_ns2(i,j,1,9) + coef_ns2(i,j,2,9) coef_ns(i,jTS,8) = coef_ns2(i,j,1,8) ! NOT added values SH j-1 coef_ns(i,jTS,7) = coef_ns2(i,j,1,7) coef_ns(i,jTS,6) = coef_ns2(i,j,1,6) coef_ns(i,jTN,1) = coef_ns(i,jTS,1) ! added values NH+ SH - just take SH added values coef_ns(i,jTN,5) = coef_ns(i,jTS,5) coef_ns(i,jTN,9) = coef_ns(i,jTS,9) coef_ns(i,jTN,2) = coef_ns(i,jTS,2) !2 am_2023.03.16 ! consider the change in j direction coef_ns(i,jTN,3) = coef_ns(i,jTS,3) !3 coef_ns(i,jTN,4) = coef_ns(i,jTS,4) !4 coef_ns(i,jTN,6) = coef_ns2(i,j,2,6) !6 03.16.2023 ! NOT added values SH j-1 coef_ns(i,jTN,7) = coef_ns2(i,j,2,7) !7 ! consider the change in j direction coef_ns(i,jTN,8) = coef_ns2(i,j,2,8) !8 coef_ns(i,jTS,10) = coef_ns2(i,j,1,10) + coef_ns2(i,j,2,10) ! added values SH + NH coef_ns(i,jTN,10) = coef_ns(i,jTS,10) !230608 begin ! do j=jlatm_JT+1,nmlat_h ! am_2023.03.15 inc. equator - loop from latm_JT to one of the equator adding both hemispheres together do j=jlatm_JT+1,nmlat_h-1 ! am_2023.03.15 inc. equator - loop from latm_JT to one of the equator adding both hemispheres together !230608 end jTS = j ! overall index SH is the same since from S-pole to equator jTN = nmlat_T1-j+1 ! overall index NH from N-pole to equator ic = 1 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ic = 5 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ic = 9 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ic = 10 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ic = 2 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ! am_2023.03.15 ic = 8 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ! ic = 3 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ! ic = 7 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ! ic = 4 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ! ic = 6 coef_ns(i,jTS,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) ! coef_ns(i,jTN,9) = coef_ns(i,jTS,9) coef_ns(i,jTN,1) = coef_ns(i,jTS,1) coef_ns(i,jTN,5) = coef_ns(i,jTS,5) coef_ns(i,jTN,10)= coef_ns(i,jTS,10) coef_ns(i,jTN,7) = coef_ns(i,jTS,7) ! 3 am_2023.03.15 changed order coef_ns(i,jTN,3) = coef_ns(i,jTS,3) ! 7 coef_ns(i,jTN,8) = coef_ns(i,jTS,8) ! 2 coef_ns(i,jTN,2) = coef_ns(i,jTS,2) ! 8 coef_ns(i,jTN,6) = coef_ns(i,jTS,6) ! 4 coef_ns(i,jTN,4) = coef_ns(i,jTS,4) ! 6 ! end do ! end lat/fieldline loop ! am_2023.03.15 inc. equator done above set equatorial values (page 14 Art's notes) ! there should be just one equator value - this remains the same !230608 begin ! jTS = nmlat_h ! overall index SH is the same since from S-pole to equator j=jTN=jTS = nmlath ! ! jTN = nmlat_T1-j+1 ! overall index NH from N-pole to equator ! coef_ns(i,jTS,2) = 0. ! am_2023.03.15 ! simplify might already be zero coef_ns(i,jTS,2:4) = 0 ! coef_ns(i,jTS,3) = 0. ! coef_ns(i,jTS,4) = 0. j = nmlat_h do ic=1,10 coef_ns(i,j,ic) = coef_ns2(i,j,1,ic)+ coef_ns2(i,j,2,ic) enddo do ic=6,8 coef_ns(i,j,ic) = .5*coef_ns(i,j,ic) enddo do ic=2,4 coef_ns(i,j,ic) = coef_ns(i,j,10-ic) enddo !230608 end ! end do ! end longitude loop ! deallocate(coef,STAT=status) if(status /= 0) write(6,*) 'deallocation of coef not succesful' deallocate(coef_ns2,STAT=status) if(status /= 0) write(6,*) 'deallocation of coef_ns2 not succesful' end subroutine set_coef_ns !----------------------------------------------------------------------------- subroutine calc_fac3 ! ! calculate FAC based on readin high latitude potential ! and add to RHS ! use fieldline_p_module,only: fieldline_p,fline_p ! implicit none ! integer :: i,j,jj,ij,isn,it,istat,jS,jN real(8), allocatable,dimension(:) :: z_ns,pot4fac_ns real(8), allocatable,dimension(:,:) :: a_ns ! if (.not.allocated(z_ns)) then allocate(z_ns(nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating z_ns') endif if (.not.allocated(a_ns)) then allocate(a_ns(nlonlat,nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating a_ns') endif if (.not.allocated(pot4fac_ns)) then allocate(pot4fac_ns(nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating pot4fac_ns') endif ! ij = 0 isn = 1 z_ns = 0 a_ns = 0 pot4fac_ns = 0 ! am 11/21/2023 putting the high latitude potential in X and then use LHS to calculate RHS FAC ! do i=1,nmlon ! loop over all longitudes j=1 ! south pole value ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,j,1)%pot_hl j=nmlat_T1 ! N-pole for each i set ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,1,2)%pot_hl ! am2023.11.21 do j=2,jlatm_JT-1 ! S-pole+1 to one southward of latm_JT both hemispheres are coupled ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,j,1)%pot_hl jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj pot4fac_ns(ij) = fline_p(i,j,2)%pot_hl end do ! end lat loop j=jlatm_JT ! SH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,j,1)%pot_hl ! am2023.05.16 jj= nmlat_T1-j+1 ! NH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+jj pot4fac_ns(ij) = fline_p(i,j,2)%pot_hl ! am_2025.08 make high latitude potential equatorward of 45 deg not symmetric ! since most high-lat. potential patterns are zero in that region I could not ! test it do j=jlatm_JT+1,nmlat_h-1 ! from latm_JT to equator - symmetric solution ij = (i-1)*nmlat_T1+j ! SH !pot4fac_ns(ij) = 0.5*(fline_p(i,j,1)%pot_hl +fline_p(i,j,2)%pot_hl) pot4fac_ns(ij) = fline_p(i,j,1)%pot_hl jj= nmlat_T1-j+1 ! NH from latm_JT to equator ij = (i-1)*nmlat_T1+jj ! note that coef_ns has the direction switched in the NH !pot4fac_ns(ij) =pot4fac_ns(ij) pot4fac_ns(ij) = fline_p(i,j,2)%pot_hl end do ! end lat loop j=nmlat_h !equator - symmetric solution ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,j,1)%pot_hl ! end do ! end longitude loop ! ! put pot4fac_ns into lhs_ns*pot4fac_ns=fac ! z_ns= 0. a_ns = lhs_ns ! just in case I am not using lhs_ns z_ns = matmul(a_ns,pot4fac_ns) ! add the FAC to the rhs rhs_ns(:) = rhs_ns(:) + z_ns(:) ! it=0 do i=1,nmlon do j=1,nmlat_T1 it = it +1 if(j.gt.nmlat_h) then write(55,'(3(i7,x),4(x,e15.8))') it,i,j,z_ns(it),pot4fac_ns(it), & fline_p(i,nmlat_T1-j+1,2)%pot_hl else write(55,'(3(i7,x),4(x,e15.8))') it,i,j,z_ns(it),pot4fac_ns(it), & fline_p(i,j,1)%pot_hl end if if(j.lt.nmlat_h) then ! SH jS = j fline_p(i,jS,1)%fac_hl = z_ns(it) if(fline_p(i,jS,1)%M3(1).ne.0) then fline_p(i,jS,1)%fac_hl = fline_p(i,jS,1)%fac_hl/fline_p(i,jS,1)%M3(1) else fline_p(i,jS,1)%fac_hl = 0. endif elseif(j.gt.nmlat_h) then ! NH jN = nmlat_T1 - j+1 fline_p(i,jN,2)%fac_hl = z_ns(it) if(fline_p(i,jN,2)%M3(1).ne.0) then fline_p(i,jN,2)%fac_hl = fline_p(i,jN,2)%fac_hl/fline_p(i,jN,2)%M3(1) else fline_p(i,jN,2)%fac_hl = 0. endif elseif(j.eq.nmlat_h) then ! equator double value jS = j jN = nmlat_T1 - j+1 fline_p(i,jS,1)%fac_hl = z_ns(it) fline_p(i,jN,2)%fac_hl = z_ns(it) if(fline_p(i,jS,isn)%M3(1).ne.0) then fline_p(i,jS,1)%fac_hl = fline_p(i,jS,1)%fac_hl/fline_p(i,jS,1)%M3(1) else fline_p(i,jS,1)%fac_hl = 0. endif if(fline_p(i,jN,2)%M3(1).ne.0) then fline_p(i,jN,2)%fac_hl = fline_p(i,jN,2)%fac_hl/fline_p(i,jN,2)%M3(1) else fline_p(i,jN,2)%fac_hl = 0. endif else write(6,*) 'subroutinecalc_fac3 stop: not allowed j value j=',j exit end if enddo enddo deallocate(z_ns,stat=istat) if (istat /= 0) call shutdown('Error deallocating z_ns') deallocate(a_ns,stat=istat) if (istat /= 0) call shutdown('Error deallocating a_ns') deallocate(pot4fac_ns,stat=istat) if (istat /= 0) call shutdown('Error deallocating pot4fac_ns') end subroutine calc_fac3 !---------------------------------------------------------------------------- subroutine fac_correct ! correct the readin FAC to make sure it is zero when integrated over the whole globe use fieldline_p_module,only: fieldline_p,fline_p use fieldline_s_module,only: fieldline_s1,fline_s1 implicit none integer :: isn,i,j,im,ip,istat,ij,jj real :: sum,sumP,sumPn,sumPs,sumn,sums,corr,corrn,corrs,facArea real, allocatable:: z_ns(:) sum=0; sumP=0; corr=0. ! correction is done globally sumn=0; sumPn=0; corrn=0. ! do correction in each hemisphere sums=0; sumPs=0; corrs=0. do isn = 1,2 ! loop over both hemisphere do i=1,nmlon ! loop over all longitudes if(i == 1) then ! wrap around in longitude im = nmlon else im = i-1 endif if(i == nmlon) then ip = 1 else ip = i+1 endif do j=2,nmlat_h ! loop over all latitudes in one hemisphere not the pole (potential set later) ! corr(i) = - ZigP(i)*abs(Jmr(i))/sinI(i) * [sum_i^N Jmr(i)*area(i)] / [sum_i^N abs(Jmr(i))*ZigP(i/sinI(i))*area(i)] ! I ignore sinI(i) by assimung it is close to 1 at high latitude ! S1 points are not really the location of P points could do an lon. average if(fline_s1(i,j,isn)%zigP.gt.0.0) then facArea = fline_p(i,j,isn)%fac_hl*fline_p(i,j,isn)%M3(1) sum = sum + facArea sumP = sumP + fline_s1(i,j,isn)%zigP*abs(facArea) ! if(isn.eq.1) then sums = sums+facArea sumPs = sumPs + fline_s1(i,j,isn)%zigP*abs(facArea) elseif(isn.eq.2) then sumn = sumn+facArea sumPn = sumPn + fline_s1(i,j,isn)%zigP*abs(facArea) endif else fline_p(i,j,isn)%fac_hl = 0. end if end do ! end lat/fieldline loop end do ! end longitude loop end do ! end hemisphere loop ! next 4 lines commented out to make sure no high latitude forcing is included ! fline_p(:,:,:)%fac_hl= 0. ! corr = 0. ! corrn= 0. ! corrs= 0. ! corr = sum/sumP ! corr(i) = [sum_i^N Jmr(i)*area(i)] / [sum_i^N abs(Jmr(i))*ZigP(i)*area(i)] corrn = sumn/sumPn ! northern hemisphere corrs = sums/sumPs ! southern hemisphere write(6,*) 'fac_hl max and min', maxval(fline_p(:,:,:)%fac_hl),minval(fline_p(:,:,:)%fac_hl) write(6,*) 'correction fac*area/[zigP*|fac|*area] glob, sh,nh ', corr, corrs, corrn write(6,*) 'sum fac*area', sum,'sums fac*area ', sums, 'sumn fac*area ',sumn ! ! correct fac_hl sum = 0.; sumn = 0. ; sums=0 do isn = 1,2 ! loop over both hemisphere do i=1,nmlon ! loop over all longitudes do j=2,nmlat_h ! loop over all latitudes in one hemisphere not the pole (potential set later) !write(66,'(4(x,e15.8))') fline_p(i,j,isn)%mlon_qd(1), & ! fline_p(i,j,isn)%mlat_qd(1),fac_hl(i,j,isn),fline_s1(i,j,isn)%zigP*abs(fac_hl(i,j,isn))*corr ! corr(i) = - ZigP(i)*abs(Jmr(i)) * [sum_i^N Jmr(i)*area(i)] / [sum_i^N abs(Jmr(i))*ZigP(i)*area(i)] ! fline_p(i,j,isn)%fac_hl = fline_p(i,j,isn)%fac_hl-fline_s1(i,j,isn)%zigP*abs(fline_p(i,j,isn)%fac_hl)*corr if(isn.eq.1) then fline_p(i,j,isn)%fac_hl = fline_p(i,j,isn)%fac_hl-fline_s1(i,j,isn)%zigP*abs(fline_p(i,j,isn)%fac_hl)*corrs sums = sums + fline_p(i,j,isn)%fac_hl ! for checking fac*area elseif(isn.eq.2) then fline_p(i,j,isn)%fac_hl = fline_p(i,j,isn)%fac_hl-fline_s1(i,j,isn)%zigP*abs(fline_p(i,j,isn)%fac_hl)*corrn sumn = sumn + fline_p(i,j,isn)%fac_hl ! for checking fac*area endif ! put in coef-array add high lat forcing to wind term ! coef_ns2(i,j,isn,10) = coef_ns2(i,j,isn,10)+fac_hl(i,j,isn) ! sum = sum + fline_p(i,j,isn)%fac_hl ! for checking fac*area ! end do ! end lat/fieldline loop end do ! end longitude loop end do ! end hemisphere loop write(6,*) 'after corr. sum fac_hl glb,sh,nh ', sum, sums,sumn ! ! add the FAC to the rhs if (.not.allocated(z_ns)) then allocate(z_ns(nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating z_ns') endif z_ns = 0 do i=1,nmlon ! loop over all longitudes j=1 ! south pole value ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) j=nmlat_T1 ! N-pole for each i set ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,1,2)%fac_hl*fline_p(i,1,2)%M3(1) do j=2,jlatm_JT-1 ! S-pole+1 to one southward of latm_JT both hemispheres are coupled ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj z_ns(ij) = fline_p(i,j,2)%fac_hl*fline_p(i,j,2)%M3(1) end do ! end lat loop j=jlatm_JT ! SH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) ! am2023.05.16 jj= nmlat_T1-j+1 ! NH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+jj z_ns(ij) = fline_p(i,j,2)%fac_hl*fline_p(i,j,2)%M3(1) do j=jlatm_JT+1,nmlat_h-1 ! from latm_JT to equator - symmetric solution ij = (i-1)*nmlat_T1+j ! SH z_ns(ij) = 0.5*(fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) +fline_p(i,j,2)%fac_hl*fline_p(i,j,2)%M3(1)) jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj ! note that coef_ns has the direction switched in the NH z_ns(ij) =z_ns(ij) end do ! end lat loop j=nmlat_h !equator - symmetric solution ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) end do ! end longitude loop rhs_ns(:) = rhs_ns(:) + z_ns(:) if (allocated(z_ns)) deallocate(z_ns,stat=istat) ! end subroutine fac_correct ! !----------------------------------------------------------------------------- subroutine const_rhs ! construct matrix lhs &rhs for solving ! set where the two hemispheres are connected. This is not in the 9-point stencil ! but needs to be done manually use fieldline_p_module,only: fieldline_p,fline_p use params_module,only: readin_amie,readin_fac ! implicit none ! integer :: i,j,jj,im,ip,ijm,it,ij,jt,status,isn,jN,jS,itcon,ii,k ! for solver integer:: info,ifail integer,parameter :: nrhmax = 1 integer :: ipiv(nlonlat) character, parameter :: trans='N' ! condition number real(8) :: colsum,anorm ,rcond,work(4*nlonlat) real(8) :: dlange integer :: info_c,iwork(nlonlat),istat character, parameter :: norm='1' ! am 1/2015 for testing real(8) :: pot_ns(nlonlat) real(8) :: z_ns(nlonlat),zpot_ns(nlonlat) real(8), allocatable,dimension(:,:) :: a_ns ! if (.not.allocated(a_ns)) then allocate(a_ns(nlonlat,nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating a_ns') endif if (.not.allocated(lhs_ns_wob)) then ! needed for test_loc_pot==.true. allocate(lhs_ns_wob(nlonlat,nlonlat),stat=status) if (istat /= 0) call shutdown('Error allocating lhs_ns_wob') endif ! z_ns = 0 lhs_ns = 0. lhs_ns_wob = 0. rhs_ns = 0. ij = 0 isn = 1 ! am 1/2015 test by putting analytical potential into stencil and compare RHS ! do i=1,nmlon ! loop over all longitudes ! if(i == 1) then ! wrap around in longitude im = nmlon else im = i-1 endif if(i == nmlon) then ip = 1 else ip = i+1 endif ! ! loop 1,nmlat_T1-1 ! S-pole to N-pole ! am 01.2023_pole add pole 2023.02 am solve two hemisphere ! j=1 ! south pole value no c6,c7,c8 values exists ij = (i-1)*nmlat_T1+j pot_ns(ij) = fline_p(i,j,1)%pot_test ! am2023.05.16 if(i.eq.1) then ! only for i=1 an equation rhs_ns(ij) = coef_ns(i,j,10) ! calculated in calc_coef - FAC not included it = (i-1)*nmlat_T1+j+1 ! c3 lhs_ns(ij,it) = coef_ns(i,j,3) ! Eq. 7.22 X_2(i)*Phi(i,2) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j ! c9 lhs_ns(ij,it) = coef_ns(i,j,9) ! Eq. 7.22 C_9*Phi^P lhs_ns_wob(ij,it) = lhs_ns(ij,it) ! 230615 begin Add bij effects art_2023.06.20 lhs_ns(ij,it) = lhs_ns(ij,it) - bijSum ! Eq.7.22 (C_9^P-beta)*Phi^P it = nmlat_T1 ! conjugate point: North pole for i=1 lhs_ns(ij,it) = bijSum ! Eq. 7.20 &7.22 beta*Phi^P* ! 230615 end do ii=2,nmlon ! from SUM_i^nmlon C3(i,1) Phi(i,2) -> Eq. 7.22 X_2(i)*Phi(i,2) it = (ii-1)*nmlat_T1+j+1 ! c3 lhs_ns(ij,it) = coef_ns(ii,j,3) lhs_ns_wob(ij,it) = lhs_ns(ij,it) end do else it = (i-1)*nmlat_T1+j ! c9 lhs_ns(ij,it) = 1. lhs_ns_wob(ij,it) = lhs_ns(ij,it) ii = 1 ! it = (ii-1)*nmlat_T1+j ! Phi(i,1) = Phi(1,1) am05.26.2025 This does not make sense lhs_ns(ij,it) = -1. lhs_ns_wob(ij,it) = lhs_ns(ij,it) endif j=nmlat_T1 ! N-pole for each i set Phi(i,Npole) = phi_pol ij = (i-1)*nmlat_T1+j pot_ns(ij) = fline_p(i,1,2)%pot_test ! am2023.05.16 rhs_ns(ij) = phi_pol it = (i-1)*nmlat_T1+j lhs_ns(ij,it) = 1. lhs_ns_wob(ij,it) = lhs_ns(ij,it) do j=2,jlatm_JT-1 ! S-pole+1 to one southward of latm_JT both hemispheres are coupled ij = (i-1)*nmlat_T1+j pot_ns(ij) = fline_p(i,j,1)%pot_test ! am2023.05.16 rhs_ns(ij) = coef_ns(i,j,10) it = (ip-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,1) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,2) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,3) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,4) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,5) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,6) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,7) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,8) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,9) -fline_p(i,j,2)%bij ! should be 1. lhs_ns_wob(ij,it) = coef_ns(i,j,9) itcon = (i-1)*nmlat_T1+nmlat_T1-j+1 ! conjugate point lhs_ns(ij,itcon)= fline_p(i,j,2)%bij ! b(i,j) Phi*(i,j) jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj ! note that coef_ns has the direction switched in the NH pot_ns(ij) = fline_p(i,j,2)%pot_test ! am2023.05.16 rhs_ns(ij) = coef_ns(i,jj,10) it = (ip-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,1) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,8) ! am2023.03.16 changed where I put the coefficients lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,7) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,6) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,5) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,4) ! am2023.03.16 lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,3) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,2) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,9) -fline_p(i,j,1)%bij lhs_ns_wob(ij,it) = coef_ns(i,jj,9) itcon = (i-1)*nmlat_T1+j ! +1 conjugate point am2013.04.17 + 1 is an error !230608 begin ! lhs_ns(ij,itcon)= lhs_ns(ij,itcon) + fline_p(i,j,1)%bij ! b(i,j) Phi*(i,j) lhs_ns(ij,itcon)= fline_p(i,j,1)%bij ! b(i,j) Phi*(i,j) !230608 end end do ! end lat loop j=jlatm_JT ! SH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+j pot_ns(ij) = fline_p(i,j,1)%pot_test ! am2023.05.16 rhs_ns(ij) = coef_ns(i,j,10) it = (ip-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,1) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,2) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,3) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,4) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,5) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,6) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,7) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,8) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,9) lhs_ns_wob(ij,it) = lhs_ns(ij,it) ! coef_ns has the direction changed- therefore ! org. coef_ns2 c6, c7, c8 at j-1 poleward -> j'+1 coef_ns c4, c3, c2 !230608 begin ! itcon = (im-1)*nmlat_T1+nmlat_T1-j+1-1 ! conjugate point C6 am2023.03.16 itcon = (im-1)*nmlat_T1+nmlat_T1-j+1+1 ! conjugate point C6 am2023.03.16 !230608 end lhs_ns(ij,itcon) = coef_ns(i,nmlat_T1-j+1,6) !230608 begin ! itcon = (i-1)*nmlat_T1+nmlat_T1-j+1-1 ! conjugate point C7 itcon = (i-1)*nmlat_T1+nmlat_T1-j+1+1 ! conjugate point C7 !230608 end lhs_ns(ij,itcon) = coef_ns(i,nmlat_T1-j+1,7) !230608 begin ! itcon = (ip-1)*nmlat_T1+nmlat_T1-j+1-1 ! conjugate point C8 itcon = (ip-1)*nmlat_T1+nmlat_T1-j+1+1 ! conjugate point C8 !230608 end lhs_ns(ij,itcon) = coef_ns(i,nmlat_T1-j+1,8) jj= nmlat_T1-j+1 ! NH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+jj ! note that coef_ns has the direction switched in the NH pot_ns(ij) = fline_p(i,j,2)%pot_test ! am2023.05.16 ! am2023.04.14 for bndry testing of test case to have symmetry in boundary condition ! rhs_ns(ij) = 0. ! am2023.04.14 for bndry testing ! it = (i-1)*nmlat_T1+jj ! lhs_ns(ij,it) = 1. ! itcon = (i-1)*nmlat_T1+j ! +1 am2023.04.17 the +1 was an error ! lhs_ns(ij,itcon) = -1. ! it = (ip-1)*nmlat_T1+jj ! lhs_ns(ij,it) = 0. ! coef_ns(i,jj,1) ! it = (ip-1)*nmlat_T1+jj+1 ! lhs_ns(ij,it) = 0. !coef_ns(i,jj,8) ! it = (i-1)*nmlat_T1+jj+1 ! lhs_ns(ij,it) = 0. !coef_ns(i,jj,7) ! it = (im-1)*nmlat_T1+jj+1 ! lhs_ns(ij,it) = 0. !coef_ns(i,jj,6) ! it = (im-1)*nmlat_T1+jj ! lhs_ns(ij,it) = 0. !coef_ns(i,jj,5) ! it = (im-1)*nmlat_T1+jj-1 ! lhs_ns(ij,it) = 0. !coef_ns(i,jj,4) ! it = (i-1)*nmlat_T1+jj-1 ! lhs_ns(ij,it) = 0. !coef_ns(i,jj,3) ! it = (ip-1)*nmlat_T1+jj-1 ! lhs_ns(ij,it) = 0. !coef_ns(i,jj,2) ! end am2023.04.14 for bndry testing !write(88,'(10(x,e15.8))') rhs_ns(ij),(lhs_ns(ij,k),k=1,nlonlat) rhs_ns(ij) = coef_ns(i,jj,10) ! do not delete original code: begin am2023.04.14 for bndry testing it = (ip-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,1) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,8) ! am2023.03.16 lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,7) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,6) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,5) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,4) ! am2023.03.16 lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,3) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,2) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,9) lhs_ns_wob(ij,it) = lhs_ns(ij,it) ! SH coef_ns has no switched direction j-1 -> j'-1 itcon = (im-1)*nmlat_T1+j-1 !+1 conjugate point C6 ! am2023.04.18 remove +1 lhs_ns(ij,itcon) = coef_ns(i,j,6) itcon = (i-1)*nmlat_T1+j-1 !+1 conjugate point C7 ! am2023.04.18 remove +1 lhs_ns(ij,itcon) = coef_ns(i,j,7) itcon = (ip-1)*nmlat_T1+j-1 !+1 conjugate point C8 ! am2023.04.18 remove +1 lhs_ns(ij,itcon) = coef_ns(i,j,8) ! end am2023.04.14 for bndry testing: original code do j=jlatm_JT+1,nmlat_h-1 ! from latm_JT to equator - symmetric solution ij = (i-1)*nmlat_T1+j ! SH pot_ns(ij) = fline_p(i,j,1)%pot_test ! am2023.05.16 rhs_ns(ij) = coef_ns(i,j,10) it = (ip-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,1) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,2) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,3) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,4) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,5) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,6) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,7) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,8) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,9) lhs_ns_wob(ij,it) = lhs_ns(ij,it) jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj ! note that coef_ns has the direction switched in the NH pot_ns(ij) = fline_p(i,j,2)%pot_test ! am2023.05.16 rhs_ns(ij) = coef_ns(i,jj,10) it = (ip-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,1) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,8) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,7) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj+1 lhs_ns(ij,it) = coef_ns(i,jj,6) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,5) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,4) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,3) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+jj-1 lhs_ns(ij,it) = coef_ns(i,jj,2) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+jj lhs_ns(ij,it) = coef_ns(i,jj,9) lhs_ns_wob(ij,it) = lhs_ns(ij,it) end do ! end lat loop j=nmlat_h !equator - symmetric solution ij = (i-1)*nmlat_T1+j pot_ns(ij) = fline_p(i,j,1)%pot_test ! am2023.05.16 rhs_ns(ij) = coef_ns(i,j,10) it = (ip-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,1) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,2) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,3) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j+1 lhs_ns(ij,it) = coef_ns(i,j,4) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,5) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (im-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,6) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,7) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (ip-1)*nmlat_T1+j-1 lhs_ns(ij,it) = coef_ns(i,j,8) lhs_ns_wob(ij,it) = lhs_ns(ij,it) it = (i-1)*nmlat_T1+j lhs_ns(ij,it) = coef_ns(i,j,9) lhs_ns_wob(ij,it) = lhs_ns(ij,it) !230607 test begin ! Write out non-zero lhs and rhs at one longitude and at latitudes around JT ! if (i.eq.nmlon/4) then ! write(77,'(3(a10,i10))') ' nmlat_T1=',nmlat_T1,' jlatm_JT=',jlatm_JT,' nlonlat=',nlonlat ! write(77,'(a36)') ' i j ij it lhs or rhs' ! jj = nmlat_T1-jlatm_JT+1 ! do j=1,nmlat_T1 ! if (j.ge.jlatm_JT-1.and.j.le.jlatm_JT+1) then ! ij = (i-1)*nmlat_T1+j ! do it=1,nlonlat ! if (lhs_ns(ij,it).ne.0.) then ! write(77,'(4i6,e15.8)') i,j,ij,it,lhs_ns(ij,it) ! endif ! enddo ! it ! write(77,'(4i6,e15.8)') i,j,ij,0,rhs_ns(ij) ! endif ! if (j.ge. jj -1.and.j.le. jj +1) then ! ij = (i-1)*nmlat_T1+j ! do it=1,nlonlat ! if (lhs_ns(ij,it).ne.0.) then ! write(77,'(4i6,e15.8)') i,j,ij,it,lhs_ns(ij,it) ! endif ! enddo ! it ! write(77,'(4i6,e15.8)') i,j,ij,0,rhs_ns(ij) ! endif ! enddo ! j ! endif !selected longitude !230607 test end ! end do ! end longitude loop ! determines FAC forcing and add to rhs if(test_loc_pot) then ! local test potential call set_local_testpot ! sets local test potential call cal_local_testFac ! calculates FAC (set dynamo forcing to zero) use matrix without bij call fac_correct_local ! correct FAC by a factor and set wind dynamo and high-lat forcing to zero else if(.not.readin_fac) then ! calculate FAC based on readin high latitude potential & add it to rhs - ! does not need correction since from the divergence of horizontal current call calc_fac3 else ! FAC is readin and needs correction call fac_correct end if ! readin_fac end if ! test_loc_pot ! ! deallocate(coef_ns,STAT=status) if(status /= 0) write(6,*) 'deallocation of coef_ns not succesful' ! ! put pot_ns into the lhs lhs_ns -> lhs_ns*pot_ns= rhs_test ! z_ns= 0. a_ns = lhs_ns ! am2023.04.17 test matrix solution ! pot_ns = rhs_ns ! z_ns = matmul( lhs_ns,pot_ns) ! it=0 ! do i=1,nmlon ! do j=1,nmlat_h ! it = it +1 ! write(55,'(2(i4,x),3(x,e15.8))') i,j,z_ns(it),rhs_ns(it),fline_p(i,j,isn)%pot_test ! enddo ! enddo ! ! do i=1,nlonlat ! loop over all longitudes ! do j=1,nlonlat ! pole to equator ! write(88,'(2(i4,x),1(x,e15.8))') i,j,lhs_ns(i,j) ! end do ! end lat/fieldline loop ! write(99,'(1(i4,x),1(x,e15.8))') i,rhs_ns(i) ! end do ! end longitude loop ! solve the matrix LHS X = RHS ! ! calculate norm anorm = 0.d0 do i=1,nlonlat colsum = 0.d0 do j=1,nlonlat ! compute norm-1 of A-> max_j SUM_i abs(a_ij) colsum = colsum + abs(lhs_ns(j,i)) end do anorm = max(anorm,colsum) end do write(6,*) 'anorm ',anorm ! ! Factorize A ! write (6,*) 'factorize A' call DGETRF(nlonlat,nlonlat,lhs_ns,nlonlat,ipiv,info) write(6,*) 'info ' , info ! if (info.EQ.0) then ! ! calculate condition number CALL DGECON( norm,nlonlat,lhs_ns,nlonlat,anorm,rcond,work,iwork,info_c) write(6,*) 'info_c ' , info_c write(6,*) 'condit=', 1/rcond ! ! Compute solution ! ! am 1/2015 test ! rhs_ns = z_ns ! it=0 ! do i=1,nmlon ! do j=1,nmlat_h ! it = it +1 ! if(j.eq.97.or.j.eq.98) rhs_ns(it) = z_ns(it) ! enddo ! enddo ! it=0 ! do i=1,nmlon ! do j=1,nmlat_h ! it = it +1 ! write(55,'(2(i4,x),3(x,e15.8))') i,j,z_ns(it),rhs_ns(it),fline_p(i,j,isn)%pot_test ! enddo ! enddo ! pot_ns = rhs_ns write (6,*) 'solve' call DGETRS(trans,nlonlat,nrhmax,lhs_ns,nlonlat,ipiv,rhs_ns,nlonlat,info) write (6,*) 'after solve',info ! z_ns = matmul(a_ns,rhs_ns) ! am2023.04.17 test ! pot_ns = 1. ! zpot_ns= matmul(a_ns,pot_ns) ! am2023.05.16 test use analystical potential ! ! ifail = 0 ! call X04CAF('General',' ',nlonlat,nrhmax,rhs_ns,nlonlat,'Solution(s)',ifail) ! it=0 do i=1,nmlon do j=1,nmlat_T1 it = it +1 if(j.gt.nmlat_h) then !if(j.eq.jlatm_JT.or.j.eq.nmlat_T1-jlatm_JT+1) then ! write(55,'(3(i7,x),1(x,e15.8))') it,i,j,rhs_ns(it) ! am2023.04.17 test !endif !write(55,'(3(i7,x),7(x,e15.8))') it,i,j,rhs_ns(it),pot_ns(it),z_ns(it),zpot_ns(it),a_ns(it,it), & ! fline_p(i,nmlat_T1-j+1,2)%mlat_m,fline_p(i,nmlat_T1-j+1,2)%mlon_m write(55,'(3(i7,x),4(x,e15.8))') it,i,j,z_ns(it),pot_ns(it), & !230607 test begin ! fline_p(i,nmlat_T1-j+1,2)%mlat_m,fline_p(i,nmlat_T1-j+1,2)%mlon_m rhs_ns(it),fline_p(i,nmlat_T1-j+1,2)%pot_test !230607 test end else !write(55,'(3(i7,x),7(x,e15.8))') it,i,j,rhs_ns(it),pot_ns(it),z_ns(it),zpot_ns(it),a_ns(it,it), & !fline_p(i,j,1)%mlat_m,fline_p(i,j,1)%mlon_m write(55,'(3(i7,x),4(x,e15.8))') it,i,j,z_ns(it),pot_ns(it), & !230607 test begin ! fline_p(i,j,1)%mlat_m,fline_p(i,j,1)%mlon_m rhs_ns(it),fline_p(i,j,1)%pot_test !230607 test end end if if(j.lt.nmlat_h) then ! SH jS = j fline_p(i,jS,1)%pot = rhs_ns(it) elseif(j.gt.nmlat_h) then ! NH jN = nmlat_T1 - j+1 fline_p(i,jN,2)%pot = rhs_ns(it) elseif(j.eq.nmlat_h) then ! equator double value jS = j jN = nmlat_T1 - j+1 fline_p(i,jS,1)%pot = rhs_ns(it) fline_p(i,jN,2)%pot = rhs_ns(it) else write(6,*) 'subroutine const_rhs stop: not allowed j value j=',j exit end if enddo enddo else write (6,*) 'The factor U is singular' end if ! deallocate(a_ns,stat=istat) ! if (istat /= 0) call shutdown('Error deallocating a_ns') end subroutine const_rhs !------------------------------------------------------------------------------- subroutine set_bij2 use fieldline_p_module,only: fieldline_p,fline_p ! use params_module,only: nmlon,nmlat_h, & ! dimension of bij ! jlatm_JT,bijSum,rho,rho_s, ! pc_uncoupled, ! switch to set bij to zero in polar cap. ! pccolatrad ! colatitude (radians) of pc bdy (nmlon) implicit none integer :: i,j,isn real :: delrho2,sumc3c7n,sumc3c7s,fac3,rho_pc ! b_mult is [|Phi|/Delta(Phi)]*(R/L)^2, where Phi is a characteristic potential value, ! Delta(Phi) is a characteristic allowed interhemispheric potential difference, ! R is Earth radius, and L is a characteristic N-S length scale for Phi. ! It is assumed that b_mult is similar for middle and auroral latitudes. real, parameter :: b_mult = 1.e6 ! original default 1.e3 do j=2,jlatm_JT ! jlatm_JT goes from pole to equator delrho2 = (rho_s(j,1) - rho_s(j-1,1))**2 do i=1,nmlon sumc3c7s = coef_ns2(i,j,1,3) + coef_ns2(i,j,1,7) sumc3c7n = coef_ns2(i,j,2,3) + coef_ns2(i,j,2,7) fline_p(i,j,1)%bij = b_mult*delrho2/(1./sumc3c7s + 1./sumc3c7n) fline_p(i,j,2)%bij = fline_p(i,j,1)%bij ! if(i.eq.10) write(6,*) 'bij_set ',i,j,fline_p(i,j,1)%bij,sumc3c7s,sumc3c7n end do enddo ! j=2,jlatm_JT ! set to zero in region where fieldlines are assumed to be equipotential ! should not be used (zero actually would mean the hemispheres are decoupled) fline_p(:,jlatm_JT+1:nmlat_h,1)%bij = 0. fline_p(:,jlatm_JT+1:nmlat_h,2)%bij = 0. ! the sum of the polar value is needed for the pole equation bijSum = 0. do i=1,nmlon bijSum = bijSum + fline_p(i,2,1)%bij end do ! bij is proportional to M3 (horizontal area of top element); ! area for j=1 is about 1/8 as large as area for j=2. We do not have to be ! precise when calculating bij, as long at it is identical at conjugate points. bijSum = bijSum/8.e0 ! do i=1,nmlon !! Note: fline_p(i,1,1)%bij and fline_p(i,1,2)%bij are not used, but are !! calculated here for completeness. Only bijSum is used. ! fline_p(i,1,1)%bij = bijSum/float(nmlon) ! fline_p(i,1,2)%bij = fline_p(i,1,1)%bij ! end do if (pc_uncoupled) then ! pccolatrad is polar-cap colatitude in radians, which for now is fixed ! but can be made variable wrt time and magnetic longitude in the future. do i=1,nmlon pccolatrad(i) = .25 ! Art set it to 0.25 -> 14deg ; 0.35 -> 20deg 0.45-> 26deg rho_pc = sin(pccolatrad(i)) ! Set bij to 0 within the polar caps, transitioning linearly ! to the full original value over a distance of about (1/3)pccolatrad. do j=1,jlatm_JT fac3 = 3.*rho(j,1)/rho_pc - 3. if (fac3.le.0.) then ! rho_pc >= rho fline_p(i,j,1)%bij = 0. fline_p(i,j,2)%bij = 0. elseif (fac3.lt.1.) then ! 4/3 rho_pc > rho fline_p(i,j,1)%bij = fac3*fline_p(i,j,1)%bij fline_p(i,j,2)%bij = fac3*fline_p(i,j,2)%bij ! the same at conjugate points since bij is the same ! If fac3>1 bij remains unmodified. endif ! if(i.eq.10) write(6,'("bij_fac ",2(x,i3),4(x,f10.7))') i,j,fline_p(i,j,1)%bij,fac3,pccolatrad(i)*rtd,(90.+ylatm(j,1)*rtd) enddo ! j enddo ! i bijSum = 0. endif ! pc_uncoupled end subroutine set_bij2 !----------------------------------------------------------------------------------------- subroutine set_local_testpot use params_module, only : pi,dtr use fieldline_p_module, only : fieldline_p,fline_p implicit none real,parameter :: & theta_max = 8.1096144*dtr , & ! determines regional extent of test potential mlat_m0 = 85.*dtr , & ! mag. latitude of center of test potential mlon_m0 = 0.*dtr ! mag. longitude of center of test potential real :: cosThetamax,costheta,fac,potConst integer :: isn,i,j cosThetamax = cos(theta_max) potConst = 1.e8 isn = 2 ! northern hemisphere do j=1,nmlat_h ! loop over latitudes (pole to equator) do i=1,nmlon ! loop over longitude costheta = sin(fline_p(i,j,isn)%mlat_m)*sin(mlat_m0)+ & cos(fline_p(i,j,isn)%mlat_m)*cos(mlat_m0)*cos(fline_p(i,j,isn)%mlon_m-mlon_m0) if(costheta.gt.cosThetamax) then fline_p(i,j,isn)%potTestLoc = potConst*(costheta-cosThetamax)**2 ! pot pos in one hemi, neg in other !write(6,'("potloc ",2(x,i4),6(x,e15.7))') j,i,fline_p(i,j,isn)%potTestLoc, & !(costheta-cosThetamax)**2,costheta,cosThetamax,fline_p(i,j,isn)%mlat_m,fline_p(i,j,isn)%mlon_m-mlon_m0 else fline_p(i,j,isn)%potTestLoc = 0. endif fline_p(i,j,1)%potTestLoc = -fline_p(i,j,isn)%potTestLoc ! southern hemisphere is negative of northern hemisphere enddo enddo end subroutine set_local_testpot !----------------------------------------------------------------------------------------- subroutine cal_local_testFac ! putting the test potential in X and then use LHS to calculate RHS FAC ! NOTE !!! same as calc_fac3 (could be merged by using either fline_p(i,j,:)%pot_hl or fline_p(i,j,:)%potTestLoc use fieldline_p_module, only : fieldline_p,fline_p implicit none integer :: i,j,jj,ij,isn,it,istat,jS,jN real(8), allocatable,dimension(:) :: z_ns,pot4fac_ns real(8), allocatable,dimension(:,:) :: a_ns if (.not.allocated(z_ns)) then allocate(z_ns(nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating z_ns') endif if (.not.allocated(a_ns)) then allocate(a_ns(nlonlat,nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating a_ns') endif if (.not.allocated(pot4fac_ns)) then allocate(pot4fac_ns(nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating pot4fac_ns') endif ij = 0 isn = 1 z_ns = 0 a_ns = 0 pot4fac_ns = 0 ! do i=1,nmlon ! loop over all longitudes j=1 ! south pole value ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,j,1)%potTestLoc j=nmlat_T1 ! N-pole for each i set ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,1,2)%potTestLoc do j=1,jlatm_JT-1 ! S-pole+1 to one southward of latm_JT both hemispheres are coupled ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,j,1)%potTestLoc jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj pot4fac_ns(ij) = fline_p(i,j,2)%potTestLoc end do ! end lat loop j=jlatm_JT ! SH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,j,1)%potTestLoc jj= nmlat_T1-j+1 ! NH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+jj pot4fac_ns(ij) = fline_p(i,j,2)%potTestLoc do j=jlatm_JT+1,nmlat_h-1 ! from latm_JT to equator - symmetric solution ij = (i-1)*nmlat_T1+j ! SH pot4fac_ns(ij) = 0.5*(fline_p(i,j,1)%potTestLoc +fline_p(i,j,2)%potTestLoc) jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj ! note that coef_ns has the direction switched in the NH pot4fac_ns(ij) =pot4fac_ns(ij) end do ! end lat loop j=nmlat_h !equator - symmetric solution ij = (i-1)*nmlat_T1+j pot4fac_ns(ij) = fline_p(i,j,1)%potTestLoc ! end do ! end longitude loop ! put pot4fac_ns into lhs_ns*pot4fac_ns=fac ! z_ns= 0. a_ns = lhs_ns_wob ! lhs_ns_wob ! this is the left hand side without the b factor z_ns = matmul(a_ns,pot4fac_ns) ! it=0 do i=1,nmlon do j=1,nmlat_T1 it = it +1 if(j.lt.nmlat_h) then ! SH jS = j fline_p(i,jS,1)%fac_hl = z_ns(it) if(fline_p(i,jS,1)%M3(1).ne.0) then fline_p(i,jS,1)%fac_hl = fline_p(i,jS,1)%fac_hl/fline_p(i,jS,1)%M3(1) else fline_p(i,jS,1)%fac_hl = 0. endif elseif(j.gt.nmlat_h) then ! NH jN = nmlat_T1 - j+1 fline_p(i,jN,2)%fac_hl = z_ns(it) if(fline_p(i,jN,2)%M3(1).ne.0) then fline_p(i,jN,2)%fac_hl = fline_p(i,jN,2)%fac_hl/fline_p(i,jN,2)%M3(1) else fline_p(i,jN,2)%fac_hl = 0. endif elseif(j.eq.nmlat_h) then ! equator double value jS = j jN = nmlat_T1 - j+1 fline_p(i,jS,1)%fac_hl = z_ns(it) fline_p(i,jN,2)%fac_hl = z_ns(it) if(fline_p(i,jS,isn)%M3(1).ne.0) then fline_p(i,jS,1)%fac_hl = fline_p(i,jS,1)%fac_hl/fline_p(i,jS,1)%M3(1) else fline_p(i,jS,1)%fac_hl = 0. endif if(fline_p(i,jN,2)%M3(1).ne.0) then fline_p(i,jN,2)%fac_hl = fline_p(i,jN,2)%fac_hl/fline_p(i,jN,2)%M3(1) else fline_p(i,jN,2)%fac_hl = 0. endif else write(6,*) 'subroutine cal_local_testFac stop: not allowed j value j=',j exit end if enddo enddo deallocate(z_ns,stat=istat) if (istat /= 0) call shutdown('Error deallocating z_ns') deallocate(a_ns,stat=istat) if (istat /= 0) call shutdown('Error deallocating a_ns') deallocate(pot4fac_ns,stat=istat) if (istat /= 0) call shutdown('Error deallocating pot4fac_ns') end subroutine cal_local_testFac ! !-------------------------------------------------------------------------------- subroutine fac_correct_local ! correct FAC from local test potential - simple integrate NH and SH and calculate a factor for adjustment ! this is different from the FAC correct of high-lat. forcing use fieldline_p_module,only: fieldline_p,fline_p implicit none integer :: isn,i,j,im,ip,istat,ij,jj real :: sumn,sums,suman,sumas,corr,facArea,signfac real, allocatable:: z_ns(:) sumn=0; sums=0; suman=0; sumas=0; corr=0. ! do correction in each hemisphere do isn = 1,2 ! loop over both hemisphere do i=1,nmlon ! loop over all longitudes do j=1,nmlat_h ! loop over all latitudes in one hemisphere not the pole (potential set later) ! corr(i) = [sum_i^N Jmr(i)*area(i) ]_sh / [sum_i^N Jmr(i)*area(i) ]_nh ! sinI(i) is not considered since fac_hl is the radial component - sinI is close to 1 at high latitude ! if(isn.eq.1) then facArea = fline_p(i,j,isn)%M3(1) if(fline_p(i,j,isn)%fac_hl.ne.0) sumas = sumas+facArea sums = sums + facArea*fline_p(i,j,isn)%fac_hl elseif(isn.eq.2) then facArea = fline_p(i,j,isn)%M3(1) if(fline_p(i,j,isn)%fac_hl.ne.0) suman = suman+facArea sumn = sumn + facArea*fline_p(i,j,isn)%fac_hl endif !if(i.eq.2) write(6,*) 'm3 ',j,fline_p(i,j,1)%M3(1),fline_p(i,j,2)%M3(1) end do ! end lat/fieldline loop end do ! end longitude loop end do ! end hemisphere loop ! sumn+corr*sums = 0 -> corr = -sumn/sums corr = -sumn/sums fline_p(:,:,1)%fac_hl = fline_p(:,:,1)%fac_hl*corr ! write(6,*) 'before correction: corr, sumn, sums, arean, areas ', corr, sumn, sums, suman, sumas ! ! check corrected fac_hl sumn = 0. ; sums = 0. do isn = 1,2 ! loop over both hemisphere do i=1,nmlon ! loop over all longitudes do j=1,nmlat_h ! loop over all latitudes in one hemisphere not the pole (potential set later) if(isn.eq.1) then facArea = fline_p(i,j,isn)%M3(1) sums = sums + facArea*fline_p(i,j,isn)%fac_hl elseif(isn.eq.2) then facArea = fline_p(i,j,isn)%M3(1) sumn = sumn + facArea*fline_p(i,j,isn)%fac_hl endif ! end do ! end lat/fieldline loop end do ! end longitude loop end do ! end hemisphere loop write(6,*) 'after correction: sumn/sums, sumn, sums ', sumn/sums, sumn, sums ! ! correct the SH test potential with the factor based on the test FAC fline_p(:,:,1)%potTestLoc = fline_p(:,:,1)%potTestLoc*corr ! ! calculate the bij*delta Phi current & subtract/add from the current at the top ! this can be used at all i,j points since bij has only non-zero values where it matters ! the bij*delta Phi current has different sign in NH & SH hemisphere which is realized by taking the ! the reverse potential difference signfac = -1. ! -1 is the correct sign ! replace rhs with test FAC if (.not.allocated(z_ns)) then allocate(z_ns(nlonlat),stat=istat) if (istat /= 0) call shutdown('Error allocating z_ns') endif z_ns = 0 do i=1,nmlon ! loop over all longitudes j=1 ! south pole value ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) ! does not need the correction since bij==0 j=nmlat_T1 ! N-pole for each i set ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,1,2)%fac_hl*fline_p(i,1,2)%M3(1) ! does not need the correction since bij==0 do j=2,jlatm_JT-1 ! S-pole+1 to one southward of latm_JT both hemispheres are coupled ! polar cap does not need the correction since bij==0 but do not check lat. index for it ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) + & signfac*fline_p(i,j,1)%bij*(fline_p(i,j,1)%potTestLoc-fline_p(i,j,2)%potTestLoc) jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj z_ns(ij) = fline_p(i,j,2)%fac_hl*fline_p(i,j,2)%M3(1)+ & signfac*fline_p(i,j,1)%bij*(fline_p(i,j,2)%potTestLoc-fline_p(i,j,1)%potTestLoc) end do ! end lat loop j=jlatm_JT ! SH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) + & signfac*fline_p(i,j,1)%bij*(fline_p(i,j,1)%potTestLoc-fline_p(i,j,2)%potTestLoc) jj= nmlat_T1-j+1 ! NH latm_JT both hemispheres are coupled at j-1 ij = (i-1)*nmlat_T1+jj z_ns(ij) = fline_p(i,j,2)%fac_hl*fline_p(i,j,2)%M3(1)+ & signfac*fline_p(i,j,1)%bij*(fline_p(i,j,2)%potTestLoc-fline_p(i,j,2)%potTestLoc) do j=jlatm_JT+1,nmlat_h-1 ! from latm_JT to equator - symmetric solution ij = (i-1)*nmlat_T1+j ! SH z_ns(ij) = 0.5*(fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) +fline_p(i,j,2)%fac_hl*fline_p(i,j,2)%M3(1)) jj= nmlat_T1-j+1 ! NH from pole to latm_JT ij = (i-1)*nmlat_T1+jj ! note that coef_ns has the direction switched in the NH z_ns(ij) =z_ns(ij) end do ! end lat loop j=nmlat_h !equator - symmetric solution ij = (i-1)*nmlat_T1+j z_ns(ij) = fline_p(i,j,1)%fac_hl*fline_p(i,j,1)%M3(1) end do ! end longitude loop ! overwrite rhs (setting wind dynamo and high-lat. fac to zero) rhs_ns(:) = z_ns(:) if (allocated(z_ns)) deallocate(z_ns,stat=istat) end subroutine fac_correct_local !-------------------------------------------------------------------------------- end module coef_module !--------------------------------------------------------------------------------------------