Project

General

Profile

tra_adv_time_iter.F90

Giovanni Aloisio, 03/06/2017 12:22 PM

 
1
!!========================================================================================================================
2
!! ***  traadv kernel extracted from the NEMO software ( http://www.nemo-ocean.eu ) ***
3
!! ***  governed by the CeCILL licence (http://www.cecill.info)  ***
4
!!
5
!!                                                       IS-ENES2 - CMCC
6
!!========================================================================================================================
7
PROGRAM tra_adv
8
#include "mpif.h"
9

    
10
   REAL, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t3sn, t3ns, t3ew, t3we
11
   REAL, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: tsn 
12
   REAL, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: pun, pvn, pwn
13
   REAL, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: mydomain, zslpx, zslpy, zwx, zwy, umask, vmask, tmask, zind
14
   REAL, ALLOCATABLE, SAVE, DIMENSION(:,:)     :: ztfreez, rnfmsk, upsmsk
15
   REAL, ALLOCATABLE, SAVE, DIMENSION(:)       :: rnfmsk_z
16
   REAL                                        :: zice, zu, z0u, zzwx, zv, z0v, zzwy, ztra, zbtr, zdt, zalpha
17
   INTEGER                                     :: nbondi, nbondj, nono, noso, noea, nowe
18
   INTEGER                                     :: jpi, jpj, jpk
19
   INTEGER                                     :: myrank, p, ierr, realcomm
20
   INTEGER                                     :: px, py, cpt
21
   CHARACTER(len=10)                           :: env
22
   INTEGER                                     :: jpreci, jprecj, nreci, nrecj, nlci, nlcj
23
   INTEGER                                     :: i, hpmid, color, a, j, ji, jj, jk, jt
24
   REAL                                        :: tstart, tend
25
   INTEGER                                     :: iter
26

    
27
   CALL get_environment_variable("PROC_X", env)
28
   READ ( env, '(i10)' ) px
29
   CALL get_environment_variable("PROC_Y", env)
30
   READ ( env, '(i10)' ) py
31
   CALL get_environment_variable("JPI", env)
32
   READ ( env, '(i10)' ) jpi
33
   CALL get_environment_variable("JPJ", env)
34
   READ ( env, '(i10)' ) jpj
35
   CALL get_environment_variable("JPK", env)
36
   READ ( env, '(i10)' ) jpk
37

    
38
   CALL get_environment_variable("CPU_PER_TASK", env)
39
   READ ( env, '(i10)' ) cpt 
40

    
41
   CALL mpi_init( ierr )
42

    
43
   CALL mpi_comm_rank( MPI_COMM_WORLD, myrank, ierr )
44
   CALL mpi_comm_size( MPI_COMM_WORLD, p, ierr )
45
      WRITE (*,*) "size=",p," PROC_X=",px, "PROC_Y", py
46
      color = MOD(myrank,cpt)
47
   CALL mpi_comm_split(MPI_COMM_WORLD, color, myrank, realcomm, ierr)
48
   IF (color == 0) THEN
49
   CALL mpi_comm_rank( realcomm, myrank, ierr )
50
   CALL mpi_comm_size( realcomm, p, ierr )
51
   
52
   IF (myrank == 0) THEN
53
      WRITE (*,*) "(PROC_X, PROC_Y)(",px,",",py,")"
54
      WRITE (*,*) "(jpi, jpj, jpk)(",jpi,",",jpj,",",jpk,")"
55
   END IF
56

    
57
   IF (p .ne. px*py) THEN
58
      WRITE (*,*) "Error! the total number of processes does not match with the values of PROC_X and PROC_Y env variables"
59
      WRITE (*,*) "size=",p," PROC_X=",px, "PROC_Y", py
60
      CALL EXIT()
61
   ENDIF
62

    
63
   hpmid = myrank + 10
64
   nowe = myrank - px
65
   noea = myrank + px
66
   noso = myrank - 1
67
   nono = myrank + 1
68
   IF (py == 1) THEN
69
      nbondi = 2
70
   ELSE IF (myrank < px) THEN
71
      nbondi = -1
72
   ELSE IF (myrank >= px*(py-1)) THEN
73
      nbondi = 1
74
   ELSE
75
      nbondi = 0
76
   END IF
77
      
78
   IF (px == 1) THEN
79
      nbondj = 2
80
   ELSE IF (MOD(myrank, px) == px-1) THEN
81
      nbondj = 1
82
   ELSE IF (MOD(myrank, px) == 0) THEN
83
      nbondj = -1
84
   ELSE
85
      nbondj = 0
86
   END IF
87
     
88

    
89
   jpreci=1
90
   jprecj=1
91
   nlci = jpi
92
   nlcj = jpj
93
   nreci = 2*jpreci
94
   nrecj = 2*jprecj
95
   ALLOCATE( mydomain (jpi,jpj,jpk))
96
   ALLOCATE( zwx (jpi,jpj,jpk))
97
   ALLOCATE( zwy (jpi,jpj,jpk))
98
   ALLOCATE( zslpx (jpi,jpj,jpk))
99
   ALLOCATE( zslpy (jpi,jpj,jpk))
100
   ALLOCATE( pun (jpi,jpj,jpk))
101
   ALLOCATE( pvn (jpi,jpj,jpk))
102
   ALLOCATE( pwn (jpi,jpj,jpk))
103
   ALLOCATE( umask (jpi,jpj,jpk))
104
   ALLOCATE( vmask (jpi,jpj,jpk))
105
   ALLOCATE( tmask (jpi,jpj,jpk))
106
   ALLOCATE( zind (jpi,jpj,jpk))
107
   ALLOCATE( ztfreez (jpi,jpj))
108
   ALLOCATE( rnfmsk (jpi,jpj))
109
   ALLOCATE( upsmsk (jpi,jpj))
110
   ALLOCATE( rnfmsk_z (jpk))
111
   ALLOCATE( tsn(jpi,jpj,jpk))
112
   ALLOCATE( t3ew(jpj,jpreci,jpk,2))
113
   ALLOCATE( t3we(jpj,jpreci,jpk,2))
114
   ALLOCATE( t3ns(jpi,jprecj,jpk,2))
115
   ALLOCATE( t3sn(jpi,jprecj,jpk,2))
116

    
117

    
118

    
119
! array initialization
120

    
121
   zwx(:,:,:)   = 0.e0
122
   zwy(:,:,:)   = 0.e0
123
   zslpx(:,:,:) = 0.e0
124
   zslpy(:,:,:) = 0.e0
125
   zind(:,:,:)  = 0.e0
126

    
127
   DO jk = 1, jpk
128
      DO jj = 1, jpj
129
          DO ji = 1, jpi
130
              tmp = .1e0 + 0.5*mod(ji+jj+jk, 2)
131
              umask(ji,jj,jk) = tmp
132
              mydomain(ji,jj,jk) = tmp
133
              pun(ji,jj,jk) = tmp
134
              pvn(ji,jj,jk) = tmp
135
              pwn(ji,jj,jk) = tmp
136
              vmask(ji,jj,jk) = tmp 
137
              tsn(ji,jj,jk) = tmp
138
              tmask(ji,jj,jk) = tmp
139
          END DO
140
      END DO
141
   END DO
142

    
143
   DO jj=1, jpj
144
      DO ji=1, jpi
145
         ztfreez(ji,jj)= .1e0 + 0.5*mod(ji+jj+jk, 2)
146
         upsmsk(ji,jj)= .1e0 + 0.5*mod(ji+jj+jk, 2)
147
         rnfmsk(ji,jj) = .1e0 + 0.5*mod(ji+jj+jk, 2)
148
      END DO
149
   END DO
150

    
151
   DO jk=1, jpk
152
      rnfmsk_z(jk)=jk
153
   END DO
154

    
155

    
156
!***********************
157
!* Start of the synphony
158
!***********************
159

    
160
   call MPI_BARRIER(MPI_COMM_WORLD, ierr)
161
   tstart = MPI_Wtime()
162

    
163
   DO ITER = 1, 100
164
      DO jk = 1, jpk
165
         DO jj = 1, jpj
166
            DO ji = 1, jpi
167
               IF( tsn(ji,jj,jk) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
168
               ELSE                                               ;   zice = 0.e0
169
               ENDIF
170
               zind(ji,jj,jk) = MAX (   &
171
                    rnfmsk(ji,jj) * rnfmsk_z(jk),      & 
172
                    upsmsk(ji,jj)               ,      &
173
                    zice                               &
174
                    &                  ) * tmask(ji,jj,jk)
175
               zind(ji,jj,jk) = 1 - zind(ji,jj,jk)
176
            END DO
177
         END DO
178
      END DO
179

    
180
      zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0
181

    
182
      DO jk = 1, jpk-1
183
         DO jj = 1, jpj-1
184
            DO ji = 1, jpi-1
185
               zwx(ji,jj,jk) = umask(ji,jj,jk) * ( mydomain(ji+1,jj,jk) - mydomain(ji,jj,jk) )
186
               zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( mydomain(ji,jj+1,jk) - mydomain(ji,jj,jk) )
187
            END DO
188
         END DO
189
      END DO
190
      CALL mpp_lnk_3d(zwx)
191
      CALL mpp_lnk_3d(zwy)
192

    
193
      zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0 
194

    
195
      DO jk = 1, jpk-1
196
         DO jj = 2, jpj
197
            DO ji = 2, jpi 
198
               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
199
                    &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
200
               zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
201
                    &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
202
            END DO
203
         END DO
204
      END DO
205

    
206
      DO jk = 1, jpk-1    
207
         DO jj = 2, jpj
208
            DO ji = 2, jpi
209
               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
210
                    &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
211
                    &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
212
               zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
213
                    &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
214
                    &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
215
            END DO
216
         END DO
217
      END DO
218

    
219
      DO jk = 1, jpk-1
220
         zdt  = 1
221
         DO jj = 2, jpj-1
222
            DO ji = 2, jpi-1
223
               z0u = SIGN( 0.5, pun(ji,jj,jk) )
224
               zalpha = 0.5 - z0u
225
               zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt
226

    
227
               zzwx = mydomain(ji+1,jj,jk) + zind(ji,jj,jk) * (zu * zslpx(ji+1,jj,jk))
228
               zzwy = mydomain(ji  ,jj,jk) + zind(ji,jj,jk) * (zu * zslpx(ji  ,jj,jk))
229

    
230
               zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
231

    
232
               z0v = SIGN( 0.5, pvn(ji,jj,jk) )
233
               zalpha = 0.5 - z0v
234
               zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt
235

    
236
               zzwx = mydomain(ji,jj+1,jk) + zind(ji,jj,jk) * (zv * zslpy(ji,jj+1,jk))
237
               zzwy = mydomain(ji,jj  ,jk) + zind(ji,jj,jk) * (zv * zslpy(ji,jj  ,jk))
238

    
239
               zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
240
            END DO
241
         END DO
242
      END DO
243

    
244
      CALL mpp_lnk_3d(zwx)
245
      CALL mpp_lnk_3d(zwy)
246

    
247
      DO jk = 1, jpk-1
248
         DO jj = 2, jpj-1     
249
            DO ji = 2, jpi-1
250
               zbtr = 1.
251
               ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
252
                    &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
253
               mydomain(ji,jj,jk) = mydomain(ji,jj,jk) + ztra
254
            END DO
255
         END DO
256
      END DO
257

    
258
      zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0
259

    
260
      DO jk = 2, jpk-1   
261
         zwx(:,:,jk) = tmask(:,:,jk) * ( mydomain(:,:,jk-1) - mydomain(:,:,jk) )
262
      END DO
263

    
264
      zslpx(:,:,1) = 0.e0
265

    
266
      DO jk = 2, jpk-1    
267
         DO jj = 1, jpj
268
            DO ji = 1, jpi
269
               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
270
                    &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
271
            END DO
272
         END DO
273
      END DO
274

    
275
      DO jk = 2, jpk-1     
276
         DO jj = 1, jpj
277
            DO ji = 1, jpi
278
               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
279
                    &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
280
                    &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
281
            END DO
282
         END DO
283
      END DO
284

    
285
      zwx(:,:, 1 ) = pwn(:,:,1) * mydomain(:,:,1)
286

    
287
      zdt  = 1
288
      zbtr = 1.
289
      DO jk = 1, jpk-1
290
         DO jj = 2, jpj-1     
291
            DO ji = 2, jpi-1
292
               z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
293
               zalpha = 0.5 + z0w
294
               zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr
295

    
296
               zzwx = mydomain(ji,jj,jk+1) + zind(ji,jj,jk) * (zw * zslpx(ji,jj,jk+1))
297
               zzwy = mydomain(ji,jj,jk  ) + zind(ji,jj,jk) * (zw * zslpx(ji,jj,jk  ))
298

    
299
               zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
300
            END DO
301
         END DO
302
      END DO
303

    
304
      zbtr = 1.
305
      DO jk = 1, jpk-1
306
         DO jj = 2, jpj-1     
307
            DO ji = 2, jpi-1
308
               ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
309
               mydomain(ji,jj,jk) =  mydomain(ji,jj,jk) + ztra
310
            END DO
311
         END DO
312
      END DO
313

    
314
   END DO
315

    
316
   call MPI_BARRIER(MPI_COMM_WORLD, ierr)
317
   tend = MPI_Wtime()
318
   write(*,*) 'minval(mydomain) avg(mydomain) maxval(mydomain):', minval(mydomain), &
319
        sum(mydomain)/size(mydomain), maxval(mydomain), 'on rank ', myrank
320
   IF (myrank == 0) THEN
321
      write(*,*) 'Elapsed time ', tend-tstart, 'on rank ', myrank
322
   END IF
323
   
324
   END IF
325
   
326
   DEALLOCATE( mydomain )
327
   DEALLOCATE( zwx )
328
   DEALLOCATE( zwy )
329
   DEALLOCATE( zslpx )
330
   DEALLOCATE( zslpy )
331
   DEALLOCATE( pun )
332
   DEALLOCATE( pvn )
333
   DEALLOCATE( pwn )
334
   DEALLOCATE( umask)
335
   DEALLOCATE( vmask)
336
   DEALLOCATE( tmask)
337
   DEALLOCATE( zind )
338
   DEALLOCATE( ztfreez )
339
   DEALLOCATE( rnfmsk)
340
   DEALLOCATE( upsmsk)
341
   DEALLOCATE( rnfmsk_z)
342
   DEALLOCATE( tsn)
343
   DEALLOCATE( t3ew)
344
   DEALLOCATE( t3we)
345
   DEALLOCATE( t3ns)
346
   DEALLOCATE( t3sn)
347

    
348
   CALL mpi_finalize(ierr)
349

    
350
CONTAINS 
351

    
352
   SUBROUTINE mpp_lnk_3d( ptab )
353
      REAL, DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab
354
      !!
355
      INTEGER  ::   ji, jj, jk, jl
356
      INTEGER  ::   imigr, iihom, ijhom
357
      INTEGER  ::   ml_req1, ml_req2, ml_err
358
      INTEGER, DIMENSION(mpi_status_size) ::   ml_stat
359
      INTEGER :: iflag
360
      INTEGER :: istatus(mpi_status_size)
361
      DOUBLE PRECISION :: t1, t2
362

    
363
      SELECT CASE ( nbondi )
364
      CASE ( -1, 0, 1 )     
365
         iihom = nlci-nreci
366
         DO jl = 1, jpreci
367
            t3ew(:,jl,:,1) = ptab(jpreci+jl,:,:)
368
            t3we(:,jl,:,1) = ptab(iihom +jl,:,:)
369
         END DO
370
      END SELECT  
371
      imigr = jpreci * jpj * jpk
372
      SELECT CASE ( nbondi ) 
373
      CASE ( -1 )
374
         CALL mpi_isend( t3we(1,1,1,1), imigr, mpi_double_precision, noea , 0, realcomm, ml_req1, iflag )
375
         CALL mpi_recv( t3ew(1,1,1,2), imigr, mpi_double_precision, noea, 0, realcomm, istatus, iflag )
376

    
377
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
378
      CASE ( 0 )
379
         CALL mpi_isend( t3ew(1,1,1,1), imigr, mpi_double_precision, nowe , 0, realcomm, ml_req1, iflag )
380
         CALL mpi_isend( t3we(1,1,1,1), imigr, mpi_double_precision, noea , 0, realcomm, ml_req2, iflag )
381
         
382
         CALL mpi_recv( t3ew(1,1,1,2), imigr, mpi_double_precision, noea, 0, realcomm, istatus, iflag )
383
         CALL mpi_recv( t3we(1,1,1,2), imigr, mpi_double_precision, nowe, 0, realcomm, istatus, iflag )
384

    
385
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
386
         CALL mpi_wait(ml_req2, ml_stat, ml_err)
387
      CASE ( 1 )
388
         CALL mpi_isend( t3ew(1,1,1,1), imigr, mpi_double_precision, nowe , 0, realcomm, ml_req1, iflag )
389
         CALL mpi_recv( t3we(1,1,1,2), imigr, mpi_double_precision, nowe, 0, realcomm, istatus, iflag )
390
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
391
      END SELECT
392

    
393
      iihom = nlci-jpreci
394
      SELECT CASE ( nbondi )
395
      CASE ( -1 )
396
         DO jl = 1, jpreci
397
            ptab(iihom+jl,:,:) = t3ew(:,jl,:,2)
398
         END DO
399
      CASE ( 0 ) 
400
         DO jl = 1, jpreci
401
            ptab(jl      ,:,:) = t3we(:,jl,:,2)
402
            ptab(iihom+jl,:,:) = t3ew(:,jl,:,2)
403
         END DO
404
      CASE ( 1 )
405
         DO jl = 1, jpreci
406
            ptab(jl      ,:,:) = t3we(:,jl,:,2)
407
         END DO
408
      END SELECT
409

    
410
      IF( nbondj /= 2 ) THEN
411
         ijhom = nlcj-nrecj
412
         DO jl = 1, jprecj
413
            t3sn(:,jl,:,1) = ptab(:,ijhom +jl,:)
414
            t3ns(:,jl,:,1) = ptab(:,jprecj+jl,:)
415
         END DO
416
      ENDIF
417
      imigr = jprecj * jpi * jpk
418
      SELECT CASE ( nbondj )     
419
      CASE ( -1 )
420
         CALL mpi_isend( t3sn(1,1,1,1), imigr, mpi_double_precision, nono , 0, realcomm, ml_req1, iflag )
421
         CALL mpi_recv( t3ns(1,1,1,2), imigr, mpi_double_precision, nono, 0, realcomm, istatus, iflag )
422

    
423
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
424
      CASE ( 0 )
425
         CALL mpi_isend( t3ns(1,1,1,1), imigr, mpi_double_precision, noso , 0, realcomm, ml_req1, iflag )
426
         CALL mpi_isend( t3sn(1,1,1,1), imigr, mpi_double_precision, nono , 0, realcomm, ml_req2, iflag )
427
         CALL mpi_recv( t3ns(1,1,1,2), imigr, mpi_double_precision, nono, 0, realcomm, istatus, iflag )
428
         CALL mpi_recv( t3sn(1,1,1,2), imigr, mpi_double_precision, noso, 0, realcomm, istatus, iflag )
429

    
430
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
431
         CALL mpi_wait(ml_req2, ml_stat, ml_err)
432
      CASE ( 1 ) 
433
         CALL mpi_isend( t3ns(1,1,1,1), imigr, mpi_double_precision, noso , 0, realcomm, ml_req1, iflag )
434
         CALL mpi_recv( t3sn(1,1,1,2), imigr, mpi_double_precision, noso, 0, realcomm, istatus, iflag )
435

    
436
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
437
      END SELECT
438
      ijhom = nlcj-jprecj
439
      SELECT CASE ( nbondj )
440
      CASE ( -1 )
441
         DO jl = 1, jprecj
442
            ptab(:,ijhom+jl,:) = t3ns(:,jl,:,2)
443
         END DO
444
      CASE ( 0 ) 
445
         DO jl = 1, jprecj
446
            ptab(:,jl      ,:) = t3sn(:,jl,:,2)
447
            ptab(:,ijhom+jl,:) = t3ns(:,jl,:,2)
448
         END DO
449
      CASE ( 1 )
450
         DO jl = 1, jprecj
451
            ptab(:,jl,:) = t3sn(:,jl,:,2)
452
         END DO
453
      END SELECT
454

    
455
   END SUBROUTINE mpp_lnk_3d
456

    
457
END PROGRAM tra_adv