sam/write_planes.f at master · tslund/sam · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
      subroutine write_planes( uu, du, u_xz, u_yz, mean, work )

c ----- The fields contained in uu array are ordered as u, v, w, [T],
c ----- p, dudz, dvdz, [dTdz], where the fields involving T are not
c ----- present for unstratified cases.

      include 'sam.h'
      include 'mpif.h'

      real   uu(Nx,Ny,nzp,Lr), mean(Nze,Lu), work(Lw)
      real   du(Nx,Ny), u_xz(Nx,nzp,2*Lu), u_yz(Ny,nzp,2*Lu)
      real   jump
      integer(kind=mpi_offset_kind) :: offset, offset_0, offset1,
     &                                 NxNze, xz_block_size,
     &                                 NyNze, yz_block_size
      integer status(mpi_status_size)

      n_frame_p = n_frame_p + 1

      n_out = 2*Lu

c ----- Write the time.

      if(l_root) then
         write(10,10) nt, time
10       format(i6,1p,e12.4)
         flush(10)
      end if

c ----- Write the means.

      if(l_root) then
         write(11,rec=n_frame_p) mean
      end if

c ----- Write the xy planes.

      do k=1, nzp
         kk = izs + k-1
         do m=1, n_xy_planes
            if( kk .eq. k_xy_plane(m) ) then
               do n=1, n_out
                  irec = (n_frame_p-1)*n_out + n
                  write(20+m,rec=irec) uu(1:Nx,1:Ny,k,n)
               end do
            end if
         end do
      end do

c ----- Gather and write the xz planes data.

      n_words_out = Nx*nzp
      NxNze = int(Nx,8)*int(Nze,8)
      xz_block_size = NxNze*int(n_out*8,8)
      offset_0 = int((n_frame_p-1),8)*xz_block_size
      offset_1 = offset_0 + int((izs-1),8)*int(Nx*8,8)

      do m=1, n_xz_planes
         j = j_xz_plane(m)

         do n=1, Lu+1
            do k=1, nzp
c         write(70+myid,'(i5,1p,4e12.4)') k, uu(40,j,k,n)
               do i=1, Nx
                  u_xz(i,k,n) = uu(i,j,k,n)
               end do
            end do
         end do

         do k=1, nzp
            L=1
            do 20 n=1, Lu

               if( n .eq. 2 ) goto 20

               L = L + 1

               call y_deriv( uu(1,1,k,n), du, work )

               do i=1, Nx
                  u_xz(i,k,Lu+L) = du(i,j)
               end do

20          continue
         end do

         do n=1, n_out
            offset = offset_1 + (n-1)*NxNze*8
c      print *, 'myid, offset = ', myid, offset
      i=10
c      do k=1, nzp
c         write(70+myid,'(i5,1p,4e12.4)') k, u_xz(i,k,n)
c      end do
            call mpi_file_write_at( fh(30+m), offset, u_xz(1,1,n),
     &              n_words_out, mpi_double_precision, status, ierr )
c            call mpi_finalize(ierr)
c            stop
         end do

      end do

c ----- Gather and write the yz planes data.

      n_words_out = Ny*nzp
      NyNze = int(Ny,8)*int(Nze,8)
      yz_block_size = NyNze*int(n_out*8,8)
      offset_0 = int((n_frame_p-1),8)*yz_block_size
      offset_1 = offset_0 + int((izs-1),8)*int(Ny*8,8)

      do m=1, n_yz_planes
         i = i_yz_plane(m)

         do n=1, Lu+1
            do k=1, nzp
               do j=1, Ny
                  u_yz(j,k,n) = uu(i,j,k,n)
               end do
            end do
         end do

         do k=1, nzp
            L=1
            do 30 n=1, Lu

               if( n .eq. 1 ) goto 30

               L = L + 1

               call x_deriv( uu(1,1,k,n), du, work )

               do j=1, Ny
                  u_yz(j,k,Lu+L) = du(i,j)
               end do

30          continue
         end do

         do n=1, n_out
            offset = offset_1 + (n-1)*NyNze*8
            call mpi_file_write_at( fh(40+m), offset, u_yz(1,1,n),
     &              n_words_out, mpi_double_precision, status, ierr )
         end do

      end do

      return
      end