[615] | 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 |
---|
[685] | 160 | SEQUENCE_C1 |
---|
[615] | 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 |
---|
[685] | 227 | |
---|
| 228 | KERNEL(remap_eta) |
---|
| 229 | ! IN : rhodz, mass_dak, mass_dbk |
---|
| 230 | ! TMP : mass_col, cur_lev, new_rhodz_cum |
---|
| 231 | ! OUT : rhodz, old_rhodz, eta |
---|
| 232 | SEQUENCE_C1 |
---|
| 233 | PROLOGUE('1') |
---|
| 234 | rhodz_cum(CELL)=0. |
---|
| 235 | cur_lev(HIDX(CELL))=1 |
---|
[687] | 236 | eta(CELL)=1. |
---|
[685] | 237 | new_rhodz_cum(HIDX(CELL))=0. |
---|
| 238 | END_BLOCK |
---|
| 239 | BODY('1,llm') |
---|
| 240 | rhodz_cum(UP(CELL)) = rhodz_cum(CELL) + rhodz(CELL) |
---|
| 241 | END_BLOCK |
---|
| 242 | EPILOGUE('llm+1') |
---|
| 243 | mass_col(HIDX(CELL)) = rhodz_cum(CELL) |
---|
| 244 | END_BLOCK |
---|
| 245 | |
---|
| 246 | BODY('1,llm') |
---|
| 247 | old_rhodz(CELL) = rhodz(CELL) |
---|
| 248 | rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL) |
---|
| 249 | rhodz_cum_target = new_rhodz_cum(HIDX(CELL)) + rhodz(CELL) |
---|
| 250 | DO level = cur_lev(HIDX(CELL)), llm |
---|
| 251 | rhodz_cum_levp1 = rhodz_cum(AT_LEVEL(CELL,level+1)) |
---|
| 252 | IF(rhodz_cum_target<=rhodz_cum_levp1) EXIT |
---|
| 253 | END DO |
---|
[687] | 254 | IF(level>llm) level=llm |
---|
[685] | 255 | rhodz_cum_lev = rhodz_cum(AT_LEVEL(CELL,level)) |
---|
| 256 | ! now rhodz_cum_lev <= rhodz_cum_target <= rhodz_cum_levp1 |
---|
| 257 | cur_lev(HIDX(CELL)) = level |
---|
| 258 | new_rhodz_cum(HIDX(CELL)) = rhodz_cum_target |
---|
| 259 | eta(UP(CELL)) = level + (rhodz_cum_target-rhodz_cum_lev)/(rhodz_cum_levp1-rhodz_cum_lev) |
---|
| 260 | END_BLOCK |
---|
| 261 | END_BLOCK |
---|
| 262 | END_BLOCK |
---|
| 263 | |
---|
| 264 | KERNEL(remap_theta) |
---|
| 265 | ! IN : thetarhodz, eta |
---|
| 266 | ! TMP : thetarhodz_cum, new_thetarhodz_cum |
---|
| 267 | ! OUT : thetarhodz |
---|
| 268 | SEQUENCE_C1 |
---|
| 269 | PROLOGUE('1') |
---|
| 270 | thetarhodz_cum(CELL)=0. |
---|
| 271 | END_BLOCK |
---|
| 272 | BODY('1,llm') |
---|
| 273 | thetarhodz_cum(UP(CELL)) = thetarhodz_cum(CELL) + thetarhodz(CELL) |
---|
| 274 | END_BLOCK |
---|
| 275 | |
---|
| 276 | BODY('1,llm+1') |
---|
| 277 | X = eta(CELL) |
---|
[687] | 278 | level = MIN(llm,FLOOR(X)) ! eta=llm+1 => level=llm, X=1 |
---|
[685] | 279 | X = X-level |
---|
| 280 | new_thetarhodz_cum(CELL) = thetarhodz_cum(AT_LEVEL(CELL,level))+X*thetarhodz(AT_LEVEL(CELL,level)) |
---|
| 281 | END_BLOCK |
---|
| 282 | BODY('1,llm') |
---|
[687] | 283 | thetarhodz(CELL) = new_thetarhodz_cum(UP(CELL)) - new_thetarhodz_cum(CELL) |
---|
[685] | 284 | END_BLOCK |
---|
| 285 | END_BLOCK |
---|
| 286 | END_BLOCK |
---|
| 287 | |
---|
| 288 | KERNEL(remap_u) |
---|
| 289 | ! IN : u, old_rhodz, rhodz, eta |
---|
| 290 | ! TMP : urhodz_cum, new_urhodz_cum |
---|
| 291 | ! OUT : u |
---|
| 292 | SEQUENCE_E0 |
---|
| 293 | PROLOGUE('1') |
---|
| 294 | urhodz_cum(EDGE)=0. |
---|
| 295 | END_BLOCK |
---|
| 296 | BODY('1,llm') |
---|
| 297 | urhodz(EDGE) = u(EDGE)*(old_rhodz(CELL1)+old_rhodz(CELL2)) |
---|
| 298 | urhodz_cum(UP(EDGE)) = urhodz_cum(EDGE) + urhodz(EDGE) |
---|
| 299 | END_BLOCK |
---|
| 300 | |
---|
| 301 | BODY('1,llm+1') |
---|
| 302 | X = .5*(eta(CELL1)+eta(CELL2)) |
---|
[687] | 303 | level = MIN(llm,FLOOR(X)) |
---|
[685] | 304 | X = X-level |
---|
| 305 | new_urhodz_cum(EDGE) = urhodz_cum(AT_LEVEL(EDGE,level))+X*urhodz(AT_LEVEL(EDGE,level)) |
---|
| 306 | END_BLOCK |
---|
| 307 | BODY('1,llm') |
---|
| 308 | u(EDGE) = (new_urhodz_cum(UP(EDGE)) - new_urhodz_cum(EDGE)) / (rhodz(CELL1)+rhodz(CELL2)) |
---|
| 309 | END_BLOCK |
---|
| 310 | END_BLOCK |
---|
| 311 | END_BLOCK |
---|
[784] | 312 | |
---|
| 313 | KERNEL(scalar_laplacian) |
---|
| 314 | FORALL_CELLS_EXT() |
---|
| 315 | ON_EDGES |
---|
| 316 | grad(EDGE) = SIGN*(b(CELL2)-b(CELL1)) |
---|
| 317 | END_BLOCK |
---|
| 318 | END_BLOCK |
---|
| 319 | |
---|
| 320 | FORALL_CELLS_EXT() |
---|
| 321 | ON_PRIMAL |
---|
| 322 | div_ij=0. |
---|
| 323 | FORALL_EDGES |
---|
| 324 | div_ij = div_ij + SIGN*LE_DE*grad(EDGE) |
---|
| 325 | END_BLOCK |
---|
| 326 | divu(CELL) = div_ij / AI |
---|
| 327 | END_BLOCK |
---|
| 328 | END_BLOCK |
---|
| 329 | |
---|
| 330 | END_BLOCK |
---|
[792] | 331 | |
---|
| 332 | KERNEL(curl_laplacian) |
---|
| 333 | FORALL_CELLS_EXT() |
---|
| 334 | ON_DUAL |
---|
| 335 | etav = 0.d0 |
---|
| 336 | FORALL_EDGES |
---|
| 337 | etav = etav + SIGN*u(EDGE) |
---|
| 338 | END_BLOCK |
---|
| 339 | qv(DUAL_CELL) = etav/AV |
---|
| 340 | END_BLOCK |
---|
| 341 | END_BLOCK |
---|
| 342 | |
---|
| 343 | FORALL_CELLS_EXT() |
---|
| 344 | ON_EDGES |
---|
| 345 | curlcurl(EDGE) = SIGN*(qv(VERTEX1)-qv(VERTEX2))*(1./LE_DE) |
---|
| 346 | END_BLOCK |
---|
| 347 | END_BLOCK |
---|
| 348 | |
---|
| 349 | END_BLOCK |
---|