1 | KERNEL(advect_vert) |
---|
2 | |
---|
3 | FORALL_CELLS_EXT('2','llm') |
---|
4 | ON_PRIMAL |
---|
5 | dzqw(CELL) = q(CELL)-q(DOWN(CELL)) |
---|
6 | END_BLOCK |
---|
7 | END_BLOCK |
---|
8 | |
---|
9 | ! We need a barrier here because we compute dzqw above and do a vertical average below |
---|
10 | BARRIER |
---|
11 | |
---|
12 | FORALL_CELLS_EXT() |
---|
13 | ON_PRIMAL |
---|
14 | CST_IFTHEN(IS_BOTTOM_LEVEL || IS_TOP_LAYER) |
---|
15 | dzq(CELL)=0. |
---|
16 | CST_ELSE |
---|
17 | dzq_down=dzqw(CELL) |
---|
18 | dzq_up=dzqw(UP(CELL) |
---|
19 | IF(dzq_up*dzq_down>0.) THEN |
---|
20 | dzqm = .5*(dzqup+dzq_down) |
---|
21 | dzqmax = max_slope * min( ABS(dzqw_down), ABS(dzq_up) ) |
---|
22 | dzq(CELL) = sign( min(abs(dzqm),dzqmax) , dzqm ) ! sign(a,b)=a*sign(b) |
---|
23 | ELSE |
---|
24 | dzq(CELL)=0. |
---|
25 | END IF |
---|
26 | CST_ENDIF |
---|
27 | END_BLOCK |
---|
28 | END_BLOCK |
---|
29 | |
---|
30 | ! We need a barrier here because we compute dzq above and do a vertical average below |
---|
31 | BARRIER |
---|
32 | |
---|
33 | FORALL_CELLS_EXT('1','llm+1') |
---|
34 | ON_PRIMAL |
---|
35 | CST_IFTHEN(IS_BOTTOM_LEVEL || IS_TOP_INTERFACE) |
---|
36 | wq(CELL)=0. |
---|
37 | CST_ELSE |
---|
38 | w = fac*wfluxt(CELL) |
---|
39 | IF(w>0.) THEN ! upward transport, upwind side is at level l |
---|
40 | sigw = w/mass(DOWN(CELL)) |
---|
41 | qq = q(DOWN(CELL) + 0.5*(1.-sigw)*dzq(DOWN(CELL)) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 |
---|
42 | ELSE ! downward transport, upwind side is at level l+1 |
---|
43 | sigw = w/mass(CELL) |
---|
44 | qq = q(CELL) -0.5*(1.+sigw)*dzq(ij,l) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0 |
---|
45 | ENDIF |
---|
46 | wq(ij,l) = w*qq |
---|
47 | CST_ENDIF |
---|
48 | END_BLOCK |
---|
49 | END_BLOCK |
---|
50 | |
---|
51 | ! We need a barrier here because we compute wq above and do a vertical difference below |
---|
52 | BARRIER |
---|
53 | |
---|
54 | FORALL_CELLS_EXT() |
---|
55 | ON_PRIMAL |
---|
56 | newmass = mass(CELL) + fac*(wflux(CELL)-wflux(UP(CELL))) |
---|
57 | q(CELL) = ( q(CELL)*mass(CELL) + wq(CELL)-wq(UP(CELL)) ) / newmass |
---|
58 | IF(update_mass) mass(CELL) = newmass |
---|
59 | END_BLOCK |
---|
60 | END_BLOCK |
---|
61 | |
---|
62 | END_BLOCK |
---|
63 | |
---|
64 | KERNEL(advect_horiz) |
---|
65 | ! evaluate tracer field at point cc using piecewise linear reconstruction |
---|
66 | ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon |
---|
67 | ! sign*hfluxt>0 iff outgoing |
---|
68 | FORALL_CELLS_EXT() |
---|
69 | ON_EDGES |
---|
70 | IF(SIGN*hfluxt(EDGE)>0.) THEN |
---|
71 | qe = qi(CELL1) |
---|
72 | qe = qe + (cc(EDGE,1)-xyz_i(HIDX(CELL1),1))*gradq3d(CELL1,1) |
---|
73 | qe = qe + (cc(EDGE,2)-xyz_i(HIDX(CELL1),2))*gradq3d(CELL1,2) |
---|
74 | qe = qe + (cc(EDGE,3)-xyz_i(HIDX(CELL1),3))*gradq3d(CELL1,3) |
---|
75 | ELSE |
---|
76 | qe = qi(CELL2) |
---|
77 | qe = qe + (cc(EDGE,1)-xyz_i(HIDX(CELL2),1))*gradq3d(CELL2,1) |
---|
78 | qe = qe + (cc(EDGE,2)-xyz_i(HIDX(CELL2),2))*gradq3d(CELL2,2) |
---|
79 | qe = qe + (cc(EDGE,3)-xyz_i(HIDX(CELL2),3))*gradq3d(CELL2,3) |
---|
80 | END IF |
---|
81 | qflux(EDGE) = hfluxt(EDGE)*qe |
---|
82 | IF(diagflux_on) qfluxt(EDGE) = qfluxt(EDGE)+qflux(EDGE) |
---|
83 | END_BLOCK |
---|
84 | END_BLOCK |
---|
85 | |
---|
86 | ! update q and, if update_mass, update mass |
---|
87 | FORALL_CELLS() |
---|
88 | ON_PRIMAL |
---|
89 | dmass=0. |
---|
90 | dq=0. |
---|
91 | FORALL_EDGES |
---|
92 | dmass = dmass + SIGN*hfluxt(EDGE) |
---|
93 | dq = dq + SIGN*qflux(EDGE) |
---|
94 | END_BLOCK |
---|
95 | newmass = mass(CELL) - dmass/AI |
---|
96 | qi(CELL) = (qi(CELL)*mass(CELL)-dq/AI) / newmass |
---|
97 | IF(update_mass) mass(CELL)=newmass |
---|
98 | END_BLOCK |
---|
99 | END_BLOCK |
---|
100 | |
---|
101 | END_BLOCK |
---|
102 | |
---|
103 | {% macro perot_reconstruction() %} |
---|
104 | ! Perot reconstruction based on Gauss theorem |
---|
105 | ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area |
---|
106 | FORALL_CELLS() |
---|
107 | ON_PRIMAL |
---|
108 | cx=centroid(HIDX(CELL),1) |
---|
109 | cy=centroid(HIDX(CELL),2) |
---|
110 | cz=centroid(HIDX(CELL),3) |
---|
111 | ux=0. ; uy=0. ; uz=0. |
---|
112 | FORALL_EDGES |
---|
113 | {{ caller() }} |
---|
114 | ux = ux + ue_le*(.5*(xyz_v(HIDX(VERTEX1),1)+xyz_v(HIDX(VERTEX2),1))-cx) |
---|
115 | uy = uy + ue_le*(.5*(xyz_v(HIDX(VERTEX1),2)+xyz_v(HIDX(VERTEX2),2))-cy) |
---|
116 | uz = uz + ue_le*(.5*(xyz_v(HIDX(VERTEX1),3)+xyz_v(HIDX(VERTEX2),3))-cz) |
---|
117 | END_BLOCK |
---|
118 | fac = scale*(1./AI) |
---|
119 | ucenter(CELL,1)=ux*fac |
---|
120 | ucenter(CELL,2)=uy*fac |
---|
121 | ucenter(CELL,3)=uz*fac |
---|
122 | END_BLOCK |
---|
123 | END_BLOCK |
---|
124 | {% endmacro %} |
---|
125 | |
---|
126 | KERNEL(wind_centered) |
---|
127 | {% call() perot_reconstruction() %} |
---|
128 | ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE)) |
---|
129 | {% endcall %} |
---|
130 | END_BLOCK |
---|
131 | |
---|
132 | KERNEL(flux_centered) |
---|
133 | ! NB : here the input data is a flux and has already the factor l_e in it |
---|
134 | ! Input data is rescaled by factor scale |
---|
135 | {% call() perot_reconstruction() %} |
---|
136 | ue_le = SIGN*ue(EDGE) |
---|
137 | {% endcall %} |
---|
138 | END_BLOCK |
---|
139 | |
---|
140 | /*------------------------------------- Energetics --------------------------------*/ |
---|
141 | |
---|
142 | {% macro energy_flux(energy) %} |
---|
143 | FORALL_CELLS_EXT() |
---|
144 | ON_PRIMAL |
---|
145 | {{ caller() }} |
---|
146 | {{ energy }}(CELL) = {{ energy }}(CELL) + frac*rhodz(CELL)*energy |
---|
147 | pk(CELL) = energy |
---|
148 | END_BLOCK |
---|
149 | END_BLOCK |
---|
150 | FORALL_CELLS_EXT() |
---|
151 | ON_EDGES |
---|
152 | {{ energy }}_flux(EDGE) = {{ energy }}_flux(EDGE) + .5*massflux(EDGE)*(pk(CELL1)+pk(CELL2)) |
---|
153 | END_BLOCK |
---|
154 | END_BLOCK |
---|
155 | {% endmacro %} |
---|
156 | |
---|
157 | KERNEL(energy_fluxes) |
---|
158 | ! First diagnose geopotential and temperature, column-wise |
---|
159 | BARRIER |
---|
160 | SEQUENCE_EXT |
---|
161 | PROLOGUE(llm) |
---|
162 | pk(CELL) = ptop + .5*g*rhodz(CELL) |
---|
163 | END_BLOCK |
---|
164 | BODY('llm-1,1,-1') |
---|
165 | pk(CELL) = pk(UP(CELL)) + (.5*g)*( rhodz(CELL)+rhodz(UP(CELL)) ) |
---|
166 | END_BLOCK |
---|
167 | END_BLOCK |
---|
168 | ! NB : at this point pressure is stored in array pk |
---|
169 | ! pk then serves as buffer to store temperature |
---|
170 | SELECT CASE(caldyn_thermo) |
---|
171 | CASE(thermo_theta) |
---|
172 | GEOPOT('temp_ik') |
---|
173 | theta_ik = theta_rhodz(CELL,1)/rhodz(CELL) |
---|
174 | theta(CELL) = theta_ik |
---|
175 | temp_ik = theta_ik*(p_ik/preff)**kappa |
---|
176 | gv = (g*Rd)*temp_ik/p_ik |
---|
177 | END_GEOPOT |
---|
178 | CASE(thermo_entropy) |
---|
179 | GEOPOT('temp_ik') |
---|
180 | theta_ik = theta_rhodz(CELL,1)/rhodz(CELL) |
---|
181 | temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp) |
---|
182 | theta(CELL) = Treff*exp(theta_ik/cpp) |
---|
183 | gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p |
---|
184 | END_GEOPOT |
---|
185 | END SELECT |
---|
186 | BARRIER |
---|
187 | ! Now accumulate energies and energy fluxes |
---|
188 | ! NB : at this point temperature is stored in array pk |
---|
189 | ! pk then serves as buffer to store energy |
---|
190 | ! enthalpy |
---|
191 | {% call energy_flux('enthalpy') %} |
---|
192 | energy = cpp*pk(CELL) |
---|
193 | {% endcall %} |
---|
194 | ! potential energy |
---|
195 | {% call energy_flux('epot') %} |
---|
196 | energy = .5*(geopot(UP(CELL))+geopot(CELL)) |
---|
197 | {% endcall %} |
---|
198 | ! theta |
---|
199 | {% call energy_flux('thetat') %} |
---|
200 | energy = theta(CELL) |
---|
201 | {% endcall %} |
---|
202 | ! kinetic energy |
---|
203 | {% call energy_flux('ekin') %} |
---|
204 | energy=0.d0 |
---|
205 | FORALL_EDGES |
---|
206 | energy = energy + le(HIDX(EDGE))*de(HIDX(EDGE))*ue(EDGE)**2 |
---|
207 | END_BLOCK |
---|
208 | energy = energy * (.25/AI) |
---|
209 | {% endcall %} |
---|
210 | ! ulon |
---|
211 | {% call energy_flux('ulon') %} |
---|
212 | cx=centroid(HIDX(CELL),1) |
---|
213 | cy=centroid(HIDX(CELL),2) |
---|
214 | cz=centroid(HIDX(CELL),3) |
---|
215 | ux=0. ; uy=0. ; uz=0. |
---|
216 | FORALL_EDGES |
---|
217 | ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE)) |
---|
218 | ux = ux + ue_le*(.5*(xyz_v(HIDX(VERTEX1),1)+xyz_v(HIDX(VERTEX2),1))-cx) |
---|
219 | uy = uy + ue_le*(.5*(xyz_v(HIDX(VERTEX1),2)+xyz_v(HIDX(VERTEX2),2))-cy) |
---|
220 | uz = uz + ue_le*(.5*(xyz_v(HIDX(VERTEX1),3)+xyz_v(HIDX(VERTEX2),3))-cz) |
---|
221 | END_BLOCK |
---|
222 | ulon_i = ux*elon_i(HIDX(CELL),1) + uy*elon_i(HIDX(CELL),2) + uz*elon_i(HIDX(CELL),3) |
---|
223 | energy = ulon_i*(1./AI) |
---|
224 | {% endcall %} |
---|
225 | |
---|
226 | END_BLOCK |
---|