source: codes/icosagcm/trunk/src/kernels/advect_horiz.k90 @ 1057

Last change on this file since 1057 was 953, checked in by adurocher, 5 years ago

trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )

File size: 3.4 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- advect_horiz ----------------------------------
3
4   !$acc data create(qflux(:,:)) present(qi(:,:), cc(:,:,:), gradq3d(:,:,:), mass(:,:), hfluxt(:,:), Ai(:), xyz_i(:,:)) async
5
6   ! evaluate tracer field at point cc using piecewise linear reconstruction
7   ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon
8   ! sign*hfluxt>0 iff outgoing
9   !$acc parallel loop collapse(2) async
10   DO l = ll_begin, ll_end
11      !DIR$ SIMD
12      DO ij=ij_begin_ext, ij_end_ext
13         IF(ne_right*hfluxt(ij+u_right,l)>0.) THEN
14            ij_tmp = ij
15         ELSE
16            ij_tmp = ij+t_right
17         END IF
18
19         qe = qi(ij_tmp,l)
20         qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
21         qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
22         qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
23
24         qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
25
26         IF(ne_lup*hfluxt(ij+u_lup,l)>0.) THEN
27            ij_tmp = ij
28         ELSE
29            ij_tmp = ij+t_lup
30         END IF
31
32         qe = qi(ij_tmp,l)
33         qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
34         qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
35         qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
36
37         qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe
38
39         IF(ne_ldown*hfluxt(ij+u_ldown,l)>0.) THEN
40            ij_tmp = ij
41         ELSE
42            ij_tmp = ij+t_ldown
43         END IF
44
45         qe = qi(ij_tmp,l)
46         qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
47         qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
48         qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
49
50         qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
51      END DO
52   END DO
53
54   IF(diagflux_on) THEN
55     !$acc parallel loop collapse(2) copy(qfluxt(:,:)) async
56     DO l = ll_begin, ll_end
57        !DIR$ SIMD
58        DO ij=ij_begin_ext, ij_end_ext
59           qfluxt(ij+u_right,l) = qfluxt(ij+u_right,l)+qflux(ij+u_right,l)
60           qfluxt(ij+u_lup,l) = qfluxt(ij+u_lup,l)+qflux(ij+u_lup,l)
61           qfluxt(ij+u_ldown,l) = qfluxt(ij+u_ldown,l)+qflux(ij+u_ldown,l)
62        END DO
63     END DO
64   END IF
65
66   ! update q and, if update_mass, update mass
67   !$acc parallel loop collapse(2) async
68   DO l = ll_begin, ll_end
69      !DIR$ SIMD
70      DO ij=ij_begin, ij_end
71         dmass=0.
72         dq=0.
73         dmass = dmass + ne_rup*hfluxt(ij+u_rup,l)
74         dq = dq + ne_rup*qflux(ij+u_rup,l)
75         dmass = dmass + ne_lup*hfluxt(ij+u_lup,l)
76         dq = dq + ne_lup*qflux(ij+u_lup,l)
77         dmass = dmass + ne_left*hfluxt(ij+u_left,l)
78         dq = dq + ne_left*qflux(ij+u_left,l)
79         dmass = dmass + ne_ldown*hfluxt(ij+u_ldown,l)
80         dq = dq + ne_ldown*qflux(ij+u_ldown,l)
81         dmass = dmass + ne_rdown*hfluxt(ij+u_rdown,l)
82         dq = dq + ne_rdown*qflux(ij+u_rdown,l)
83         dmass = dmass + ne_right*hfluxt(ij+u_right,l)
84         dq = dq + ne_right*qflux(ij+u_right,l)
85         newmass = mass(ij,l) - dmass/Ai(ij)
86         qi(ij,l) = (qi(ij,l)*mass(ij,l)-dq/Ai(ij)) / newmass
87         IF(update_mass) mass(ij,l)=newmass
88      END DO
89   END DO
90   
91   !$acc end data
92
93   !---------------------------- advect_horiz ----------------------------------
94   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.