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