Project

General

Profile

tra_adv.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

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

    
36
   CALL get_environment_variable("CPU_PER_TASK", env)
37
   READ ( env, '(i10)' ) cpt 
38

    
39
   CALL mpi_init( ierr )
40

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

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

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

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

    
115

    
116

    
117
! array initialization
118

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

    
125
   DO jk = 1, jpk
126
      DO jj = 1, jpj
127
          DO ji = 1, jpi
128
              umask(ji,jj,jk) = ji*jj*jk
129
              mydomain(ji,jj,jk) =ji*jj*jk
130
              pun(ji,jj,jk) =ji*jj*jk
131
              pvn(ji,jj,jk) =ji*jj*jk
132
              pwn(ji,jj,jk) =ji*jj*jk
133
              vmask(ji,jj,jk)= ji*jj*jk
134
              tsn(ji,jj,jk)= ji*jj*jk
135
              tmask(ji,jj,jk)= ji*jj*jk
136
          END DO
137
      END DO
138
   END DO
139

    
140
   DO jj=1, jpj
141
      DO ji=1, jpi
142
         ztfreez(ji,jj)=ji*jj
143
         upsmsk(ji,jj)=ji*jj
144
         rnfmsk(ji,jj) = ji*jj
145
      END DO
146
   END DO
147

    
148
   DO jk=1, jpk
149
      rnfmsk_z(jk)=jk
150
   END DO
151

    
152

    
153
!***********************
154
!* Start of the synphony
155
!***********************
156
   DO jk = 1, jpk
157
      DO jj = 1, jpj
158
         DO ji = 1, jpi
159
            IF( tsn(ji,jj,jk) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
160
            ELSE                                               ;   zice = 0.e0
161
            ENDIF
162
            zind(ji,jj,jk) = MAX (   &
163
               rnfmsk(ji,jj) * rnfmsk_z(jk),      & 
164
               upsmsk(ji,jj)               ,      &
165
               zice                               &
166
               &                  ) * tmask(ji,jj,jk)
167
               zind(ji,jj,jk) = 1 - zind(ji,jj,jk)
168
         END DO
169
      END DO
170
   END DO
171

    
172
   zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0
173

    
174
   DO jk = 1, jpk-1
175
      DO jj = 1, jpj-1
176
         DO ji = 1, jpi-1
177
             zwx(ji,jj,jk) = umask(ji,jj,jk) * ( mydomain(ji+1,jj,jk) - mydomain(ji,jj,jk) )
178
             zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( mydomain(ji,jj+1,jk) - mydomain(ji,jj,jk) )
179
         END DO
180
      END DO
181
   END DO
182
   CALL mpp_lnk_3d(zwx)
183
   CALL mpp_lnk_3d(zwy)
184

    
185
   zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0 
186

    
187
   DO jk = 1, jpk-1
188
     DO jj = 2, jpj
189
        DO ji = 2, jpi 
190
           zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
191
           &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
192
           zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
193
           &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
194
        END DO
195
     END DO
196
  END DO
197

    
198
  DO jk = 1, jpk-1    
199
     DO jj = 2, jpj
200
        DO ji = 2, jpi
201
           zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
202
           &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
203
           &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
204
           zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
205
           &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
206
           &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
207
        END DO
208
     END DO
209
  END DO 
210

    
211
  DO jk = 1, jpk-1
212
     zdt  = 1
213
     DO jj = 2, jpj-1
214
        DO ji = 2, jpi-1
215
            z0u = SIGN( 0.5, pun(ji,jj,jk) )
216
            zalpha = 0.5 - z0u
217
            zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt
218

    
219
            zzwx = mydomain(ji+1,jj,jk) + zind(ji,jj,jk) * (zu * zslpx(ji+1,jj,jk))
220
            zzwy = mydomain(ji  ,jj,jk) + zind(ji,jj,jk) * (zu * zslpx(ji  ,jj,jk))
221

    
222
            zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
223
            
224
            z0v = SIGN( 0.5, pvn(ji,jj,jk) )
225
            zalpha = 0.5 - z0v
226
            zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt
227

    
228
            zzwx = mydomain(ji,jj+1,jk) + zind(ji,jj,jk) * (zv * zslpy(ji,jj+1,jk))
229
            zzwy = mydomain(ji,jj  ,jk) + zind(ji,jj,jk) * (zv * zslpy(ji,jj  ,jk))
230

    
231
            zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
232
         END DO
233
      END DO
234
  END DO
235

    
236
  CALL mpp_lnk_3d(zwx)
237
  CALL mpp_lnk_3d(zwy)
238

    
239
  DO jk = 1, jpk-1
240
     DO jj = 2, jpj-1     
241
        DO ji = 2, jpi-1
242
           zbtr = 1.
243
           ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
244
           &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
245
           mydomain(ji,jj,jk) = mydomain(ji,jj,jk) + ztra
246
        END DO
247
     END DO
248
  END DO
249

    
250
  zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0
251

    
252
  DO jk = 2, jpk-1   
253
     zwx(:,:,jk) = tmask(:,:,jk) * ( mydomain(:,:,jk-1) - mydomain(:,:,jk) )
254
  END DO
255

    
256
  zslpx(:,:,1) = 0.e0
257

    
258
  DO jk = 2, jpk-1    
259
     DO jj = 1, jpj
260
        DO ji = 1, jpi
261
           zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
262
           &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
263
        END DO
264
     END DO
265
  END DO
266

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

    
277
  zwx(:,:, 1 ) = pwn(:,:,1) * mydomain(:,:,1)
278

    
279
  zdt  = 1
280
  zbtr = 1.
281
  DO jk = 1, jpk-1
282
     DO jj = 2, jpj-1     
283
        DO ji = 2, jpi-1
284
           z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
285
           zalpha = 0.5 + z0w
286
           zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr
287

    
288
           zzwx = mydomain(ji,jj,jk+1) + zind(ji,jj,jk) * (zw * zslpx(ji,jj,jk+1))
289
           zzwy = mydomain(ji,jj,jk  ) + zind(ji,jj,jk) * (zw * zslpx(ji,jj,jk  ))
290

    
291
           zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
292
        END DO
293
     END DO
294
  END DO
295

    
296
  zbtr = 1.
297
  DO jk = 1, jpk-1
298
     DO jj = 2, jpj-1     
299
        DO ji = 2, jpi-1
300
           ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
301
           mydomain(ji,jj,jk) =  mydomain(ji,jj,jk) + ztra
302
        END DO
303
     END DO
304
  END DO
305

    
306
  END IF
307

    
308
  DEALLOCATE( mydomain )
309
  DEALLOCATE( zwx )
310
  DEALLOCATE( zwy )
311
  DEALLOCATE( zslpx )
312
  DEALLOCATE( zslpy )
313
  DEALLOCATE( pun )
314
  DEALLOCATE( pvn )
315
  DEALLOCATE( pwn )
316
  DEALLOCATE( umask)
317
  DEALLOCATE( vmask)
318
  DEALLOCATE( tmask)
319
  DEALLOCATE( zind )
320
  DEALLOCATE( ztfreez )
321
  DEALLOCATE( rnfmsk)
322
  DEALLOCATE( upsmsk)
323
  DEALLOCATE( rnfmsk_z)
324
  DEALLOCATE( tsn)
325
  DEALLOCATE( t3ew)
326
  DEALLOCATE( t3we)
327
  DEALLOCATE( t3ns)
328
  DEALLOCATE( t3sn)
329

    
330
   CALL mpi_finalize(ierr)
331

    
332
CONTAINS 
333

    
334
   SUBROUTINE mpp_lnk_3d( ptab )
335
      REAL, DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab
336
      !!
337
      INTEGER  ::   ji, jj, jk, jl
338
      INTEGER  ::   imigr, iihom, ijhom
339
      INTEGER  ::   ml_req1, ml_req2, ml_err
340
      INTEGER, DIMENSION(mpi_status_size) ::   ml_stat
341
      INTEGER :: iflag
342
      INTEGER :: istatus(mpi_status_size)
343
      DOUBLE PRECISION :: t1, t2
344

    
345
      SELECT CASE ( nbondi )
346
      CASE ( -1, 0, 1 )     
347
         iihom = nlci-nreci
348
         DO jl = 1, jpreci
349
            t3ew(:,jl,:,1) = ptab(jpreci+jl,:,:)
350
            t3we(:,jl,:,1) = ptab(iihom +jl,:,:)
351
         END DO
352
      END SELECT  
353
      imigr = jpreci * jpj * jpk
354
      SELECT CASE ( nbondi ) 
355
      CASE ( -1 )
356
         CALL mpi_isend( t3we(1,1,1,1), imigr, mpi_double_precision, noea , 0, realcomm, ml_req1, iflag )
357
         CALL mpi_recv( t3ew(1,1,1,2), imigr, mpi_double_precision, noea, 0, realcomm, istatus, iflag )
358

    
359
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
360
      CASE ( 0 )
361
         CALL mpi_isend( t3ew(1,1,1,1), imigr, mpi_double_precision, nowe , 0, realcomm, ml_req1, iflag )
362
         CALL mpi_isend( t3we(1,1,1,1), imigr, mpi_double_precision, noea , 0, realcomm, ml_req2, iflag )
363
         
364
         CALL mpi_recv( t3ew(1,1,1,2), imigr, mpi_double_precision, noea, 0, realcomm, istatus, iflag )
365
         CALL mpi_recv( t3we(1,1,1,2), imigr, mpi_double_precision, nowe, 0, realcomm, istatus, iflag )
366

    
367
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
368
         CALL mpi_wait(ml_req2, ml_stat, ml_err)
369
      CASE ( 1 )
370
         CALL mpi_isend( t3ew(1,1,1,1), imigr, mpi_double_precision, nowe , 0, realcomm, ml_req1, iflag )
371
         CALL mpi_recv( t3we(1,1,1,2), imigr, mpi_double_precision, nowe, 0, realcomm, istatus, iflag )
372
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
373
      END SELECT
374

    
375
      iihom = nlci-jpreci
376
      SELECT CASE ( nbondi )
377
      CASE ( -1 )
378
         DO jl = 1, jpreci
379
            ptab(iihom+jl,:,:) = t3ew(:,jl,:,2)
380
         END DO
381
      CASE ( 0 ) 
382
         DO jl = 1, jpreci
383
            ptab(jl      ,:,:) = t3we(:,jl,:,2)
384
            ptab(iihom+jl,:,:) = t3ew(:,jl,:,2)
385
         END DO
386
      CASE ( 1 )
387
         DO jl = 1, jpreci
388
            ptab(jl      ,:,:) = t3we(:,jl,:,2)
389
         END DO
390
      END SELECT
391

    
392
      IF( nbondj /= 2 ) THEN
393
         ijhom = nlcj-nrecj
394
         DO jl = 1, jprecj
395
            t3sn(:,jl,:,1) = ptab(:,ijhom +jl,:)
396
            t3ns(:,jl,:,1) = ptab(:,jprecj+jl,:)
397
         END DO
398
      ENDIF
399
      imigr = jprecj * jpi * jpk
400
      SELECT CASE ( nbondj )     
401
      CASE ( -1 )
402
         CALL mpi_isend( t3sn(1,1,1,1), imigr, mpi_double_precision, nono , 0, realcomm, ml_req1, iflag )
403
         CALL mpi_recv( t3ns(1,1,1,2), imigr, mpi_double_precision, nono, 0, realcomm, istatus, iflag )
404

    
405
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
406
      CASE ( 0 )
407
         CALL mpi_isend( t3ns(1,1,1,1), imigr, mpi_double_precision, noso , 0, realcomm, ml_req1, iflag )
408
         CALL mpi_isend( t3sn(1,1,1,1), imigr, mpi_double_precision, nono , 0, realcomm, ml_req2, iflag )
409
         CALL mpi_recv( t3ns(1,1,1,2), imigr, mpi_double_precision, nono, 0, realcomm, istatus, iflag )
410
         CALL mpi_recv( t3sn(1,1,1,2), imigr, mpi_double_precision, noso, 0, realcomm, istatus, iflag )
411

    
412
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
413
         CALL mpi_wait(ml_req2, ml_stat, ml_err)
414
      CASE ( 1 ) 
415
         CALL mpi_isend( t3ns(1,1,1,1), imigr, mpi_double_precision, noso , 0, realcomm, ml_req1, iflag )
416
         CALL mpi_recv( t3sn(1,1,1,2), imigr, mpi_double_precision, noso, 0, realcomm, istatus, iflag )
417

    
418
         CALL mpi_wait(ml_req1, ml_stat, ml_err)
419
      END SELECT
420
      ijhom = nlcj-jprecj
421
      SELECT CASE ( nbondj )
422
      CASE ( -1 )
423
         DO jl = 1, jprecj
424
            ptab(:,ijhom+jl,:) = t3ns(:,jl,:,2)
425
         END DO
426
      CASE ( 0 ) 
427
         DO jl = 1, jprecj
428
            ptab(:,jl      ,:) = t3sn(:,jl,:,2)
429
            ptab(:,ijhom+jl,:) = t3ns(:,jl,:,2)
430
         END DO
431
      CASE ( 1 )
432
         DO jl = 1, jprecj
433
            ptab(:,jl,:) = t3sn(:,jl,:,2)
434
         END DO
435
      END SELECT
436

    
437
   END SUBROUTINE mpp_lnk_3d
438

    
439
END PROGRAM tra_adv