source: codes/icosagcm/devel/Python/src/kernels_transport.jin @ 792

Last change on this file since 792 was 792, checked in by jisesh, 6 years ago

devel/unstructured : added kernel for curl curl ; used by Baroclinic_3D_ullrich
File size: 10.0 KB
Line 
1KERNEL(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
62END_BLOCK
63
64KERNEL(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
101END_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
126KERNEL(wind_centered)
127{% call() perot_reconstruction() %}
128ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE))
129{% endcall %}
130END_BLOCK
131
132KERNEL(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() %}
136ue_le = SIGN*ue(EDGE)
137{% endcall %}
138END_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
157KERNEL(energy_fluxes)
158  ! First diagnose geopotential and temperature, column-wise
159  BARRIER
160  SEQUENCE_C1
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
226END_BLOCK
227
228KERNEL(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
236        eta(CELL)=1.
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
254        IF(level>llm) level=llm
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
262END_BLOCK
263
264KERNEL(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)
278        level = MIN(llm,FLOOR(X)) ! eta=llm+1 => level=llm, X=1
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')
283        thetarhodz(CELL) = new_thetarhodz_cum(UP(CELL)) - new_thetarhodz_cum(CELL)
284      END_BLOCK
285    END_BLOCK
286END_BLOCK
287
288KERNEL(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))
303        level = MIN(llm,FLOOR(X))
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
311END_BLOCK
312
313KERNEL(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
330END_BLOCK
331
332KERNEL(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
349END_BLOCK
Note: See TracBrowser for help on using the repository browser.