[850] | 1 | MODULE compute_caldyn_slow_NH_mod |
---|
[939] | 2 | USE grid_param |
---|
[850] | 3 | IMPLICIT NONE |
---|
| 4 | PRIVATE |
---|
| 5 | |
---|
[876] | 6 | #include "../unstructured/unstructured.h90" |
---|
| 7 | |
---|
[850] | 8 | PUBLIC :: compute_caldyn_slow_NH |
---|
| 9 | |
---|
| 10 | CONTAINS |
---|
| 11 | |
---|
[876] | 12 | #ifdef BEGIN_DYSL |
---|
| 13 | |
---|
| 14 | KERNEL(caldyn_slow_NH) |
---|
| 15 | FORALL_CELLS_EXT('1', 'llm+1') |
---|
| 16 | ON_PRIMAL |
---|
| 17 | CST_IF(IS_INNER_INTERFACE, w_il(CELL) = 2.*W(CELL)/(rhodz(KDOWN(CELL))+rhodz(KUP(CELL))) ) |
---|
| 18 | CST_IF(IS_BOTTOM_LEVEL, w_il(CELL) = 2.*W(CELL)/rhodz(KUP(CELL)) ) |
---|
| 19 | CST_IF(IS_TOP_INTERFACE, w_il(CELL) = 2.*W(CELL)/rhodz(KDOWN(CELL)) ) |
---|
| 20 | END_BLOCK |
---|
| 21 | END_BLOCK |
---|
| 22 | FORALL_CELLS_EXT('1', 'llm+1') |
---|
| 23 | ON_EDGES |
---|
| 24 | ! compute DePhi, v_el, G_el, F_el |
---|
| 25 | ! v_el, W2_el and therefore G_el incorporate metric factor le_de |
---|
| 26 | ! while DePhil, W_el and F_el do not |
---|
| 27 | W_el = .5*( W(CELL2)+W(CELL1) ) |
---|
| 28 | DePhil(EDGE) = SIGN*(Phi(CELL2)-Phi(CELL1)) |
---|
| 29 | F_el(EDGE) = DePhil(EDGE)*W_el |
---|
| 30 | W2_el = .5*LE_DE * ( W(CELL1)*w_il(CELL1) + W(CELL2)*w_il(CELL2) ) |
---|
| 31 | v_el(EDGE) = .5*LE_DE*(u(KUP(EDGE))+u(KDOWN(EDGE))) ! checked |
---|
| 32 | G_el(EDGE) = v_el(EDGE)*W_el - DePhil(EDGE)*W2_el |
---|
| 33 | END_BLOCK |
---|
| 34 | END_BLOCK |
---|
| 35 | |
---|
| 36 | FORALL_CELLS_EXT('1', 'llm+1') |
---|
| 37 | ! compute GradPhi2, dPhi, dW |
---|
| 38 | ON_PRIMAL |
---|
| 39 | gPhi2=0. |
---|
| 40 | dP=0. |
---|
| 41 | divG=0 |
---|
| 42 | FORALL_EDGES |
---|
| 43 | gPhi2 = gPhi2 + LE_DE*DePhil(EDGE)**2 |
---|
| 44 | dP = dP + LE_DE*DePhil(EDGE)*v_el(EDGE) |
---|
| 45 | divG = divG + SIGN*G_el(EDGE) ! -div(G_el), G_el already has le_de |
---|
| 46 | END_BLOCK |
---|
| 47 | gradPhi2(CELL) = 1./(2.*AI) * gPhi2 |
---|
| 48 | dPhi(CELL) = gradPhi2(CELL)*w_il(CELL) - 1./(2.*AI)*dP |
---|
| 49 | dW(CELL) = (-1./AI)*divG |
---|
| 50 | END_BLOCK |
---|
| 51 | END_BLOCK |
---|
| 52 | |
---|
| 53 | ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below |
---|
| 54 | BARRIER |
---|
| 55 | |
---|
| 56 | FORALL_CELLS_EXT() |
---|
| 57 | ! Compute berni at scalar points |
---|
| 58 | ON_PRIMAL |
---|
| 59 | u2=0. |
---|
| 60 | FORALL_EDGES |
---|
| 61 | u2 = u2 + LE_DE*u(EDGE)**2 |
---|
| 62 | END_BLOCK |
---|
| 63 | berni(CELL) = 1./(4.*AI) * u2 - .25*( gradPhi2(CELL)*w_il(CELL)**2 + gradPhi2(UP(CELL))*w_il(UP(CELL))**2 ) |
---|
| 64 | END_BLOCK |
---|
| 65 | END_BLOCK |
---|
| 66 | |
---|
| 67 | FORALL_CELLS_EXT() |
---|
| 68 | ON_EDGES |
---|
| 69 | ! Compute mass flux and grad(berni) |
---|
| 70 | uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE) - .5*( F_el(EDGE)+F_el(UP(EDGE)) ) |
---|
| 71 | hflux(EDGE) = LE_DE*uu |
---|
| 72 | du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2)) |
---|
| 73 | END_BLOCK |
---|
| 74 | END_BLOCK |
---|
| 75 | |
---|
| 76 | END_BLOCK |
---|
| 77 | |
---|
| 78 | #endif END_DYSL |
---|
| 79 | |
---|
| 80 | SUBROUTINE compute_caldyn_vert_NH_unst(mass,geopot,W,wflux, eta_dot,wcov,W_etadot, du,dPhi,dW) |
---|
| 81 | USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT |
---|
[939] | 82 | USE data_unstructured_mod, ONLY : enter_trace, exit_trace, & |
---|
| 83 | id_vert_NH, left,right |
---|
[876] | 84 | FIELD_MASS :: mass, eta_dot, wcov, W_etadot ! IN, BUF*3 |
---|
| 85 | FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW ! IN*3, INOUT*2 |
---|
| 86 | FIELD_U :: du ! INOUT |
---|
| 87 | DECLARE_INDICES |
---|
| 88 | NUM :: w_ij, wflux_ij |
---|
| 89 | START_TRACE(id_vert_NH, 6,0,1) |
---|
| 90 | #include "../kernels_unst/caldyn_vert_NH.k90" |
---|
| 91 | STOP_TRACE |
---|
| 92 | END SUBROUTINE compute_caldyn_vert_NH_unst |
---|
| 93 | |
---|
[850] | 94 | SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW) |
---|
| 95 | USE icosa |
---|
[851] | 96 | USE trace |
---|
[850] | 97 | USE caldyn_vars_mod |
---|
| 98 | USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1 |
---|
| 99 | REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) ! prognostic "velocity" |
---|
| 100 | REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) ! rho*dz |
---|
| 101 | REAL(rstd),INTENT(IN) :: Phi(iim*jjm,llm+1) ! prognostic geopotential |
---|
| 102 | REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1) ! prognostic vertical momentum |
---|
| 103 | |
---|
| 104 | REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s |
---|
| 105 | REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm) |
---|
| 106 | REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1) |
---|
| 107 | REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1) |
---|
| 108 | |
---|
| 109 | REAL(rstd) :: w_il(iim*jjm,llm+1) ! Wil/mil |
---|
| 110 | REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux |
---|
| 111 | REAL(rstd) :: gradPhi2(iim*jjm,llm+1) ! grad_Phi**2 |
---|
| 112 | REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi) |
---|
| 113 | |
---|
| 114 | INTEGER :: ij,l,kdown,kup |
---|
| 115 | REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu |
---|
| 116 | |
---|
| 117 | REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function |
---|
| 118 | REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W |
---|
| 119 | REAL(rstd) :: v_el(3*iim*jjm,llm+1) |
---|
| 120 | |
---|
| 121 | REAL(rstd) :: berni1(iim*jjm) ! Bernoulli function |
---|
| 122 | REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W |
---|
| 123 | REAL(rstd) :: v_el1(3*iim*jjm) |
---|
| 124 | |
---|
| 125 | CALL trace_start("compute_caldyn_slow_NH") |
---|
| 126 | |
---|
| 127 | IF(dysl) THEN |
---|
| 128 | |
---|
| 129 | !$OMP BARRIER |
---|
| 130 | #include "../kernels_hex/caldyn_slow_NH.k90" |
---|
| 131 | !$OMP BARRIER |
---|
| 132 | |
---|
| 133 | ELSE |
---|
| 134 | |
---|
| 135 | #define BERNI(ij) berni1(ij) |
---|
| 136 | #define G_EL(ij) G_el1(ij) |
---|
| 137 | #define V_EL(ij) v_el1(ij) |
---|
| 138 | |
---|
| 139 | DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces) |
---|
| 140 | IF(l==1) THEN |
---|
| 141 | kdown=1 |
---|
| 142 | ELSE |
---|
| 143 | kdown=l-1 |
---|
| 144 | END IF |
---|
| 145 | IF(l==llm+1) THEN |
---|
| 146 | kup=llm |
---|
| 147 | ELSE |
---|
| 148 | kup=l |
---|
| 149 | END IF |
---|
| 150 | ! below : "checked" means "formula also valid when kup=kdown (top/bottom)" |
---|
| 151 | ! compute mil, wil=Wil/mil |
---|
| 152 | DO ij=ij_begin_ext, ij_end_ext |
---|
| 153 | w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked |
---|
| 154 | END DO |
---|
| 155 | ! compute DePhi, v_el, G_el, F_el |
---|
| 156 | ! v_el, W2_el and therefore G_el incorporate metric factor le_de |
---|
| 157 | ! while DePhil, W_el and F_el don't |
---|
| 158 | DO ij=ij_begin_ext, ij_end_ext |
---|
| 159 | ! Compute on edge 'right' |
---|
| 160 | W_el = .5*( W(ij,l)+W(ij+t_right,l) ) |
---|
| 161 | DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l)) |
---|
| 162 | F_el(ij+u_right,l) = DePhil(ij+u_right,l)*W_el |
---|
| 163 | W2_el = .5*le_de(ij+u_right) * & |
---|
| 164 | ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) ) |
---|
| 165 | V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked |
---|
| 166 | G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el |
---|
| 167 | ! Compute on edge 'lup' |
---|
| 168 | W_el = .5*( W(ij,l)+W(ij+t_lup,l) ) |
---|
| 169 | DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l)) |
---|
| 170 | F_el(ij+u_lup,l) = DePhil(ij+u_lup,l)*W_el |
---|
| 171 | W2_el = .5*le_de(ij+u_lup) * & |
---|
| 172 | ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) ) |
---|
| 173 | V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked |
---|
| 174 | G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el |
---|
| 175 | ! Compute on edge 'ldown' |
---|
| 176 | W_el = .5*( W(ij,l)+W(ij+t_ldown,l) ) |
---|
| 177 | DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l)) |
---|
| 178 | F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el |
---|
| 179 | W2_el = .5*le_de(ij+u_ldown) * & |
---|
| 180 | ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) ) |
---|
| 181 | V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked |
---|
| 182 | G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el |
---|
| 183 | END DO |
---|
| 184 | ! compute GradPhi2, dPhi, dW |
---|
| 185 | DO ij=ij_begin_ext, ij_end_ext |
---|
| 186 | gradPhi2(ij,l) = & |
---|
| 187 | 1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + & |
---|
| 188 | le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 + & |
---|
| 189 | le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 + & |
---|
| 190 | le_de(ij+u_left)*DePhil(ij+u_left,l)**2 + & |
---|
| 191 | le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 + & |
---|
| 192 | le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 ) |
---|
| 193 | |
---|
| 194 | dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & |
---|
| 195 | ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi, |
---|
| 196 | DePhil(ij+u_rup,l)*V_EL(ij+u_rup) + & ! v_el already has le_de |
---|
| 197 | DePhil(ij+u_lup,l)*V_EL(ij+u_lup) + & |
---|
| 198 | DePhil(ij+u_left,l)*V_EL(ij+u_left) + & |
---|
| 199 | DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + & |
---|
| 200 | DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) ) |
---|
| 201 | |
---|
| 202 | dW(ij,l) = -1./Ai(ij)*( & ! -div(G_el), |
---|
| 203 | ne_right*G_EL(ij+u_right) + & ! G_el already has le_de |
---|
| 204 | ne_rup*G_EL(ij+u_rup) + & |
---|
| 205 | ne_lup*G_EL(ij+u_lup) + & |
---|
| 206 | ne_left*G_EL(ij+u_left) + & |
---|
| 207 | ne_ldown*G_EL(ij+u_ldown) + & |
---|
| 208 | ne_rdown*G_EL(ij+u_rdown)) |
---|
| 209 | END DO |
---|
| 210 | END DO |
---|
| 211 | |
---|
| 212 | DO l=ll_begin, ll_end ! compute on k levels (layers) |
---|
| 213 | ! Compute berni at scalar points |
---|
| 214 | DO ij=ij_begin_ext, ij_end_ext |
---|
| 215 | BERNI(ij) = & |
---|
| 216 | 1/(4*Ai(ij))*( & |
---|
| 217 | le_de(ij+u_right)*u(ij+u_right,l)**2 + & |
---|
| 218 | le_de(ij+u_rup)*u(ij+u_rup,l)**2 + & |
---|
| 219 | le_de(ij+u_lup)*u(ij+u_lup,l)**2 + & |
---|
| 220 | le_de(ij+u_left)*u(ij+u_left,l)**2 + & |
---|
| 221 | le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & |
---|
| 222 | le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) & |
---|
| 223 | - .25*( gradPhi2(ij,l) *w_il(ij,l)**2 + & |
---|
| 224 | gradPhi2(ij,l+1)*w_il(ij,l+1)**2 ) |
---|
| 225 | END DO |
---|
| 226 | ! Compute mass flux and grad(berni) at edges |
---|
| 227 | DO ij=ij_begin_ext, ij_end_ext |
---|
| 228 | ! Compute on edge 'right' |
---|
| 229 | uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) & |
---|
| 230 | -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) |
---|
| 231 | hflux(ij+u_right,l) = uu_right*le_de(ij+u_right) |
---|
| 232 | du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right)) |
---|
| 233 | ! Compute on edge 'lup' |
---|
| 234 | uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) & |
---|
| 235 | -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1)) |
---|
| 236 | hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup) |
---|
| 237 | du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup)) |
---|
| 238 | ! Compute on edge 'ldown' |
---|
| 239 | uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) & |
---|
| 240 | -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) |
---|
| 241 | hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown) |
---|
| 242 | du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) |
---|
| 243 | END DO |
---|
| 244 | END DO |
---|
| 245 | |
---|
| 246 | #undef V_EL |
---|
| 247 | #undef G_EL |
---|
| 248 | #undef BERNI |
---|
| 249 | |
---|
| 250 | END IF ! dysl |
---|
| 251 | |
---|
| 252 | CALL trace_end("compute_caldyn_slow_NH") |
---|
| 253 | |
---|
| 254 | END SUBROUTINE compute_caldyn_slow_NH |
---|
| 255 | |
---|
| 256 | END MODULE compute_caldyn_slow_NH_mod |
---|