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

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

trunk : backported commits r582-r598 (transport diagnostics)
File size: 3.7 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(diagflux_on) qfluxt(ij+u_right,l) = qfluxt(ij+u_right,l)+qflux(ij+u_right,l)
22         IF(ne_lup*hfluxt(ij+u_lup,l)>0.) THEN
23            qe = qi(ij,l)
24            qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
25            qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
26            qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
27         ELSE
28            qe = qi(ij+t_lup,l)
29            qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij+t_lup,1))*gradq3d(ij+t_lup,l,1)
30            qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij+t_lup,2))*gradq3d(ij+t_lup,l,2)
31            qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij+t_lup,3))*gradq3d(ij+t_lup,l,3)
32         END IF
33         qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe
34         IF(diagflux_on) qfluxt(ij+u_lup,l) = qfluxt(ij+u_lup,l)+qflux(ij+u_lup,l)
35         IF(ne_ldown*hfluxt(ij+u_ldown,l)>0.) THEN
36            qe = qi(ij,l)
37            qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
38            qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
39            qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
40         ELSE
41            qe = qi(ij+t_ldown,l)
42            qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij+t_ldown,1))*gradq3d(ij+t_ldown,l,1)
43            qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij+t_ldown,2))*gradq3d(ij+t_ldown,l,2)
44            qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij+t_ldown,3))*gradq3d(ij+t_ldown,l,3)
45         END IF
46         qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
47         IF(diagflux_on) qfluxt(ij+u_ldown,l) = qfluxt(ij+u_ldown,l)+qflux(ij+u_ldown,l)
48      END DO
49   END DO
50   ! update q and, if update_mass, update mass
51   DO l = ll_begin, ll_end
52      !DIR$ SIMD
53      DO ij=ij_begin, ij_end
54         dmass=0.
55         dq=0.
56         dmass = dmass + ne_rup*hfluxt(ij+u_rup,l)
57         dq = dq + ne_rup*qflux(ij+u_rup,l)
58         dmass = dmass + ne_lup*hfluxt(ij+u_lup,l)
59         dq = dq + ne_lup*qflux(ij+u_lup,l)
60         dmass = dmass + ne_left*hfluxt(ij+u_left,l)
61         dq = dq + ne_left*qflux(ij+u_left,l)
62         dmass = dmass + ne_ldown*hfluxt(ij+u_ldown,l)
63         dq = dq + ne_ldown*qflux(ij+u_ldown,l)
64         dmass = dmass + ne_rdown*hfluxt(ij+u_rdown,l)
65         dq = dq + ne_rdown*qflux(ij+u_rdown,l)
66         dmass = dmass + ne_right*hfluxt(ij+u_right,l)
67         dq = dq + ne_right*qflux(ij+u_right,l)
68         newmass = mass(ij,l) - dmass/Ai(ij)
69         qi(ij,l) = (qi(ij,l)*mass(ij,l)-dq/Ai(ij)) / newmass
70         IF(update_mass) mass(ij,l)=newmass
71      END DO
72   END DO
73   !---------------------------- advect_horiz ----------------------------------
74   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.