surprise

src/inital.F90

@@ -66,7 +66,13 @@ SUBROUTINE inital
       CALL update_ghosts_moments
       CALL solve_EM_fields
       CALL update_ghosts_EM
-    END SELECT
+    CASE('ricci')
+      CALL speak('Init Ricci')
+      CALL init_ricci ! init only gyrocenter density
+      CALL update_ghosts_moments
+      CALL solve_EM_fields
+      CALL update_ghosts_EM
+  END SELECT
   ENDIF
   ! closure of j>J, p>P and j<0, p<0 moments
   CALL speak('Apply closure')
@@ -448,3 +454,83 @@ SUBROUTINE init_ppj
   ENDDO
 END SUBROUTINE init_ppj
 !******************************************************************************!
+
+!******************************************************************************!
+!!!!!!! Initialize Ricci density
+!******************************************************************************!
+SUBROUTINE init_ricci
+  USE grid,       ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
+                        local_nkx_offset, local_nky_offset,&
+                        ngp, ngj, ngz, iky0, parray, jarray, contains_ky0, AA_x, AA_y
+  USE fields,     ONLY: moments
+  USE prec_const, ONLY: xp, imagu
+  USE initial_par,ONLY: iseed, init_noiselvl, init_background
+  USE model,      ONLY: LINEARITY
+  USE parallel,   ONLY: my_id
+  IMPLICIT NONE
+
+  COMPLEX(xp), DIMENSION(186,52) :: ricci_mat_real, ricci_mat_imag, ricci_face
+  REAL(xp) :: tmp_real, tmp_imag
+  COMPLEX(xp) :: scaling
+  INTEGER  :: ia,ip,ij,ikx,iky,iz
+  ! open data file
+  ricci_mat_real = 0; ricci_mat_imag = 0
+  open(unit = 1 , file = "/home/ahoffman/gyacomo/Gallery/fourier_ricci_real.txt")
+  read(1,*) ricci_mat_real
+  close(1)
+  open(unit = 2 , file = "/home/ahoffman/gyacomo/Gallery/fourier_ricci_imag.txt")
+  read(2,*) ricci_mat_imag
+  close(2)
+  scaling = 0.000002
+  moments = 0._xp
+    !**** Broad noise initialization *******************************************
+  DO ia=1,local_na
+    DO ikx=1,total_nkx
+      DO iky=1,local_nky
+        DO ip=1+ngp/2,local_np+ngp/2
+          DO ij=1+ngj/2,local_nj+ngj/2
+            DO iz=1+ngz/2,local_nz+ngz/2
+              IF((ikx+local_nkx_offset .LE. 186) .AND. (iky+local_nky_offset .LE. 52)) THEN
+                IF ( (parray(ip) .EQ. 0) .AND. (jarray(ij) .EQ. 0) ) THEN
+                  moments(ia,ip,ij,iky,ikx,iz,:) = scaling*(ricci_mat_real(ikx+local_nkx_offset,iky+local_nky_offset)&
+                  + imagu*ricci_mat_imag(ikx+local_nkx_offset,iky+local_nky_offset))
+                ELSE
+                  moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
+                ENDIF
+              ELSE
+                moments(ia,ip,ij,iky,ikx,iz,:) = 0._xp
+              ENDIF
+            END DO
+          END DO
+        END DO
+      END DO
+    END DO
+    print*, sum(moments)
+    IF ( contains_ky0 ) THEN
+      DO ip=1+ngp/2,local_np+ngp/2
+        DO ij=1+ngj/2,local_nj+ngj/2
+          DO iz=1+ngz/2,local_nz+ngz/2
+            DO ikx=2,total_nkx/2 !symmetry at ky = 0 for all z
+              moments(ia, ip,ij,iky0,ikx,:,:) = moments(ia, ip,ij,iky0,total_nkx+2-ikx,:,:)
+            END DO
+          END DO
+        END DO
+      END DO
+    ENDIF
+    ! Putting to zero modes that are not in the 2/3 Orszag rule
+    IF (LINEARITY .NE. 'linear') THEN
+      DO ikx=1,total_nkx
+        DO iky=1,local_nky
+          DO iz=1,local_nz+ngz
+            DO ip=1,local_np+ngp
+              DO ij=1,local_nj+ngj
+                moments(ia, ip,ij,iky,ikx,iz, :) = moments(ia, ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
+              ENDDO
+            ENDDO
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDIF
+  ENDDO
+END SUBROUTINE init_ricci
+!******************************************************************************!
\ No newline at end of file