1 | KERNEL(caldyn_vert) |
---|
2 | SEQUENCE |
---|
3 | BODY('llm-1,1,-1') |
---|
4 | ! cumulate mass flux convergence from top to bottom |
---|
5 | convm(CELL) = convm(CELL) + convm(UP(CELL)) |
---|
6 | END_BLOCK |
---|
7 | EPILOGUE(1) |
---|
8 | dmass_col(HIDX(CELL)) = convm(CELL) |
---|
9 | END_BLOCK |
---|
10 | BODY('2,llm') |
---|
11 | ! Compute vertical mass flux (l=1,llm+1 set to zero at init) |
---|
12 | wflux(CELL) = mass_bl(CELL) * dmass_col(HIDX(CELL)) - convm(CELL) |
---|
13 | END_BLOCK |
---|
14 | END_BLOCK |
---|
15 | |
---|
16 | ! make sure wflux is up to date |
---|
17 | BARRIER |
---|
18 | |
---|
19 | FORALL_CELLS('2','llm') |
---|
20 | ON_EDGES |
---|
21 | wwuu(EDGE) = .5*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)-u(DOWN(EDGE))) |
---|
22 | END_BLOCK |
---|
23 | END_BLOCK |
---|
24 | |
---|
25 | ! make sure wwuu is up to date |
---|
26 | BARRIER |
---|
27 | |
---|
28 | FORALL_CELLS() |
---|
29 | ON_PRIMAL |
---|
30 | convm(CELL) = mass_dbk(CELL) * dmass_col(HIDX(CELL)) ! FIXME : we shoud prognose mass_col |
---|
31 | END_BLOCK |
---|
32 | ON_EDGES |
---|
33 | du(EDGE) = du(EDGE) - (wwuu(EDGE)+wwuu(UP(EDGE))) / (rhodz(CELL1)+rhodz(CELL2)) |
---|
34 | END_BLOCK |
---|
35 | END_BLOCK |
---|
36 | |
---|
37 | DO iq=1,nqdyn |
---|
38 | FORALL_CELLS('2', 'llm') |
---|
39 | ON_PRIMAL |
---|
40 | dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) + 0.5*(theta(CELL,iq)+theta(DOWN(CELL),iq))*wflux(CELL) |
---|
41 | END_BLOCK |
---|
42 | END_BLOCK |
---|
43 | FORALL_CELLS('1', 'llm-1') |
---|
44 | ON_PRIMAL |
---|
45 | dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) - 0.5*(theta(CELL,iq)+theta(UP(CELL),iq))*wflux(UP(CELL)) |
---|
46 | END_BLOCK |
---|
47 | END_BLOCK |
---|
48 | END DO |
---|
49 | END_BLOCK |
---|
50 | |
---|
51 | KERNEL(gradient) |
---|
52 | FORALL_CELLS_EXT() |
---|
53 | ON_EDGES |
---|
54 | grad(EDGE) = SIGN*(b(CELL2)-b(CELL1)) |
---|
55 | END_BLOCK |
---|
56 | END_BLOCK |
---|
57 | END_BLOCK |
---|
58 | |
---|
59 | KERNEL(div) |
---|
60 | FORALL_CELLS_EXT() |
---|
61 | ON_PRIMAL |
---|
62 | div_ij=0. |
---|
63 | FORALL_EDGES |
---|
64 | div_ij = div_ij + SIGN*LE_DE*u(EDGE) |
---|
65 | END_BLOCK |
---|
66 | divu(CELL) = div_ij / AI |
---|
67 | END_BLOCK |
---|
68 | END_BLOCK |
---|
69 | END_BLOCK |
---|
70 | |
---|
71 | KERNEL(theta) |
---|
72 | IF(caldyn_eta==eta_mass) THEN ! Compute mass |
---|
73 | ! FIXME : here mass_col is computed from rhodz |
---|
74 | ! so that the DOFs are the same whatever caldyn_eta |
---|
75 | ! in DYNAMICO mass_col is prognosed rather than rhodz |
---|
76 | SEQUENCE_EXT |
---|
77 | PROLOGUE(0) |
---|
78 | mass_col(HIDX(CELL))=0. |
---|
79 | END_BLOCK |
---|
80 | BODY('1,llm') |
---|
81 | mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL) |
---|
82 | END_BLOCK |
---|
83 | END_BLOCK |
---|
84 | FORALL_CELLS_EXT() |
---|
85 | ON_PRIMAL |
---|
86 | ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass |
---|
87 | ! m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL) |
---|
88 | ! rhodz(CELL) = m/g |
---|
89 | rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL) |
---|
90 | END_BLOCK |
---|
91 | END_BLOCK |
---|
92 | END IF |
---|
93 | DO iq=1,nqdyn |
---|
94 | FORALL_CELLS_EXT() |
---|
95 | ON_PRIMAL |
---|
96 | theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL) |
---|
97 | END_BLOCK |
---|
98 | END_BLOCK |
---|
99 | END DO |
---|
100 | END_BLOCK |
---|
101 | |
---|
102 | KERNEL(caldyn_fast) |
---|
103 | ! |
---|
104 | SELECT CASE(caldyn_thermo) |
---|
105 | CASE(thermo_boussinesq) |
---|
106 | FORALL_CELLS() |
---|
107 | ON_PRIMAL |
---|
108 | berni(CELL) = pk(CELL) |
---|
109 | ! from now on pk contains the vertically-averaged geopotential |
---|
110 | pk(CELL) = .5*(geopot(CELL)+geopot(UP(CELL))) |
---|
111 | END_BLOCK |
---|
112 | END_BLOCK |
---|
113 | CASE(thermo_theta) |
---|
114 | FORALL_CELLS() |
---|
115 | ON_PRIMAL |
---|
116 | berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL))) |
---|
117 | END_BLOCK |
---|
118 | END_BLOCK |
---|
119 | CASE(thermo_entropy) |
---|
120 | FORALL_CELLS() |
---|
121 | ON_PRIMAL |
---|
122 | berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL))) |
---|
123 | berni(CELL) = berni(CELL) + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s) |
---|
124 | END_BLOCK |
---|
125 | END_BLOCK |
---|
126 | CASE DEFAULT |
---|
127 | PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo ! FIXME |
---|
128 | STOP |
---|
129 | END SELECT |
---|
130 | ! |
---|
131 | FORALL_CELLS() |
---|
132 | ON_EDGES |
---|
133 | due = .5*(theta(CELL1,1)+theta(CELL2,1))*(pk(CELL2)-pk(CELL1)) + berni(CELL2)-berni(CELL1) |
---|
134 | du(EDGE) = du(EDGE) - SIGN*due |
---|
135 | u(EDGE) = u(EDGE) + tau*du(EDGE) |
---|
136 | END_BLOCK |
---|
137 | END_BLOCK |
---|
138 | ! |
---|
139 | END_BLOCK |
---|
140 | |
---|
141 | KERNEL(pvort_only) |
---|
142 | FORALL_CELLS_EXT() |
---|
143 | ON_DUAL |
---|
144 | etav = 0.d0 |
---|
145 | FORALL_EDGES |
---|
146 | etav = etav + SIGN*u(EDGE) |
---|
147 | END_BLOCK |
---|
148 | hv=0. |
---|
149 | FORALL_VERTICES |
---|
150 | hv = hv + RIV2*rhodz(VERTEX) |
---|
151 | END_BLOCK |
---|
152 | qv(DUAL_CELL) = (etav + FV*AV )/(hv*AV) |
---|
153 | END_BLOCK |
---|
154 | END_BLOCK |
---|
155 | |
---|
156 | FORALL_CELLS() |
---|
157 | ON_EDGES |
---|
158 | qu(EDGE)=0.5d0*(qv(VERTEX1)+qv(VERTEX2)) |
---|
159 | END_BLOCK |
---|
160 | END_BLOCK |
---|
161 | END_BLOCK |
---|
162 | |
---|
163 | KERNEL(coriolis) |
---|
164 | ! |
---|
165 | DO iq=1,nqdyn |
---|
166 | FORALL_CELLS_EXT() |
---|
167 | ON_EDGES |
---|
168 | Ftheta(EDGE) = .5*(theta(CELL1,iq)+theta(CELL2,iq))*hflux(EDGE) |
---|
169 | END_BLOCK |
---|
170 | END_BLOCK |
---|
171 | FORALL_CELLS() |
---|
172 | ON_PRIMAL |
---|
173 | divF=0. |
---|
174 | FORALL_EDGES |
---|
175 | divF = divF + Ftheta(EDGE)*SIGN |
---|
176 | END_BLOCK |
---|
177 | dtheta_rhodz(CELL,iq) = -divF / AI |
---|
178 | END_BLOCK |
---|
179 | END_BLOCK |
---|
180 | END DO ! iq |
---|
181 | ! |
---|
182 | FORALL_CELLS() |
---|
183 | ON_PRIMAL |
---|
184 | divF=0. |
---|
185 | FORALL_EDGES |
---|
186 | divF = divF + hflux(EDGE)*SIGN |
---|
187 | END_BLOCK |
---|
188 | convm(CELL) = -divF / AI |
---|
189 | END_BLOCK |
---|
190 | END_BLOCK |
---|
191 | ! |
---|
192 | FORALL_CELLS() |
---|
193 | ON_EDGES |
---|
194 | du_trisk=0. |
---|
195 | FORALL_TRISK |
---|
196 | du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK)) |
---|
197 | END_BLOCK |
---|
198 | du(EDGE) = du(EDGE) + .5*du_trisk |
---|
199 | END_BLOCK |
---|
200 | END_BLOCK |
---|
201 | |
---|
202 | END_BLOCK |
---|