source: codes/icosagcm/devel/src/kernels/advect_horiz.k90 @ 582

Last change on this file since 582 was 582, checked in by dubos, 7 years ago

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