[853] | 1 | MODULE compute_caldyn_solver_mod |
---|
| 2 | USE grid_param, ONLY : llm |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | PRIVATE |
---|
| 5 | |
---|
[878] | 6 | #include "../unstructured/unstructured.h90" |
---|
| 7 | |
---|
[853] | 8 | PUBLIC :: compute_caldyn_solver |
---|
| 9 | |
---|
| 10 | CONTAINS |
---|
| 11 | |
---|
[878] | 12 | #ifdef BEGIN_DYSL |
---|
| 13 | |
---|
| 14 | KERNEL(caldyn_solver) |
---|
| 15 | ! |
---|
| 16 | ! Compute pressure (pres) and Exner function (pk) |
---|
| 17 | ! kappa = R/Cp |
---|
| 18 | ! 1-kappa = Cv/Cp |
---|
| 19 | ! Cp/Cv = 1/(1-kappa) |
---|
| 20 | gamma = 1./(1.-kappa) |
---|
| 21 | vreff = Rd*Treff/preff ! reference specific volume |
---|
| 22 | Cvd = 1./(cpp-Rd) |
---|
| 23 | Rd_preff = kappa*cpp/preff |
---|
| 24 | FORALL_CELLS_EXT() |
---|
| 25 | SELECT CASE(caldyn_thermo) |
---|
| 26 | CASE(thermo_theta) |
---|
| 27 | ON_PRIMAL |
---|
| 28 | rho_ij = 1./(geopot(UP(CELL))-geopot(CELL)) |
---|
| 29 | rho_ij = rho_ij*g*rhodz(CELL) |
---|
| 30 | X_ij = Rd_preff*theta(CELL,1)*rho_ij |
---|
| 31 | ! kappa.theta.rho = p/exner |
---|
| 32 | ! => X = (p/p0)/(exner/Cp) |
---|
| 33 | ! = (p/p0)^(1-kappa) |
---|
| 34 | pres(CELL) = preff*(X_ij**gamma) ! pressure |
---|
| 35 | ! Compute Exner function (needed by compute_caldyn_fast) |
---|
| 36 | ! other formulae possible if exponentiation is slow |
---|
| 37 | pk(CELL) = cpp*((pres(CELL)/preff)**kappa) ! Exner |
---|
| 38 | END_BLOCK |
---|
| 39 | CASE(thermo_entropy) |
---|
| 40 | ON_PRIMAL |
---|
| 41 | rho_ij = 1./(geopot(UP(CELL))-geopot(CELL)) |
---|
| 42 | rho_ij = rho_ij*g*rhodz(CELL) |
---|
| 43 | T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd ) |
---|
| 44 | pres(CELL) = rho_ij*Rd*T_ij |
---|
| 45 | pk(CELL) = T_ij |
---|
| 46 | END_BLOCK |
---|
| 47 | CASE DEFAULT |
---|
| 48 | STOP |
---|
| 49 | END SELECT |
---|
| 50 | END_BLOCK |
---|
| 51 | |
---|
| 52 | ! We need a barrier here because we compute pres above and do a vertical difference below |
---|
| 53 | BARRIER |
---|
| 54 | |
---|
| 55 | FORALL_CELLS_EXT('1', 'llm+1') |
---|
| 56 | ON_PRIMAL |
---|
| 57 | CST_IFTHEN(IS_BOTTOM_LEVEL) |
---|
| 58 | ! Lower BC |
---|
| 59 | dW(CELL) = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL) |
---|
| 60 | CST_ELSEIF(IS_TOP_INTERFACE) |
---|
| 61 | ! Top BC |
---|
| 62 | dW(CELL) = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL) |
---|
| 63 | CST_ELSE |
---|
| 64 | dW(CELL) = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL) |
---|
| 65 | CST_ENDIF |
---|
| 66 | W(CELL) = W(CELL)+tau*dW(CELL) ! update W |
---|
| 67 | dPhi(CELL) = g*g*W(CELL)/m_il(CELL) |
---|
| 68 | END_BLOCK |
---|
| 69 | END_BLOCK |
---|
| 70 | |
---|
| 71 | ! We need a barrier here because we update W above and do a vertical average below |
---|
| 72 | BARRIER |
---|
| 73 | |
---|
| 74 | FORALL_CELLS_EXT() |
---|
| 75 | ON_PRIMAL |
---|
| 76 | ! compute du = -0.5*g^2.grad(w^2) |
---|
| 77 | berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 ) |
---|
| 78 | END_BLOCK |
---|
| 79 | END_BLOCK |
---|
| 80 | FORALL_CELLS_EXT() |
---|
| 81 | ON_EDGES |
---|
| 82 | du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2)) |
---|
| 83 | END_BLOCK |
---|
| 84 | END_BLOCK |
---|
| 85 | END_BLOCK |
---|
| 86 | |
---|
| 87 | #endif END_DYSL |
---|
| 88 | |
---|
| 89 | SUBROUTINE compute_caldyn_solver_unst(tau,rhodz,theta, berni,pres,m_il, pk,geopot,W,dPhi,dW,du) |
---|
| 90 | USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT |
---|
| 91 | USE earth_const |
---|
| 92 | USE trace |
---|
| 93 | USE grid_param, ONLY : nqdyn |
---|
| 94 | USE disvert_mod, ONLY : ptop |
---|
| 95 | USE data_unstructured_mod, ONLY : id_solver,primal_num,dual_num,edge_num,left, right,PHI_BOT, & |
---|
| 96 | enter_trace, exit_trace |
---|
| 97 | USE compute_NH_geopot_mod, ONLY : compute_NH_geopot_unst |
---|
| 98 | NUM, INTENT(IN) :: tau |
---|
| 99 | FIELD_MASS :: rhodz,pk,berni,pres ! IN, OUT, BUF*2 |
---|
| 100 | FIELD_THETA :: theta ! IN |
---|
| 101 | FIELD_GEOPOT :: geopot,W,dPhi,dW, m_il ! INOUT,INOUT, OUT,OUT, BUF |
---|
| 102 | FIELD_U :: du ! OUT |
---|
| 103 | DECLARE_INDICES |
---|
| 104 | NUM :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff, Rd_preff |
---|
| 105 | REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6 ! FIXME |
---|
| 106 | #define PHI_BOT(ij) Phi_bot |
---|
| 107 | #include "../kernels_unst/caldyn_mil.k90" |
---|
| 108 | IF(tau>0) THEN ! solve implicit problem for geopotential |
---|
| 109 | CALL compute_NH_geopot_unst(tau, rhodz, m_il, theta, W, geopot) |
---|
| 110 | END IF |
---|
| 111 | START_TRACE(id_solver, 7,0,1) |
---|
| 112 | #include "../kernels_unst/caldyn_solver.k90" |
---|
| 113 | STOP_TRACE |
---|
| 114 | #undef PHI_BOT |
---|
| 115 | END SUBROUTINE compute_caldyn_solver_unst |
---|
| 116 | |
---|
[853] | 117 | SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du) |
---|
| 118 | USE icosa |
---|
| 119 | USE caldyn_vars_mod |
---|
| 120 | USE trace |
---|
| 121 | USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1 |
---|
| 122 | USE disvert_mod, ONLY : ptop |
---|
| 123 | USE caldyn_kernels_hevi_mod |
---|
[859] | 124 | USE compute_NH_geopot_mod |
---|
[853] | 125 | REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6 |
---|
| 126 | REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs |
---|
| 127 | REAL(rstd),INTENT(IN) :: phis(iim*jjm) |
---|
| 128 | REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) |
---|
| 129 | REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn) |
---|
| 130 | REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm) |
---|
| 131 | REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) |
---|
| 132 | REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0 |
---|
| 133 | REAL(rstd),INTENT(OUT) :: m_il(iim*jjm,llm+1) ! rhodz averaged to interfaces |
---|
| 134 | REAL(rstd),INTENT(OUT) :: pres(iim*jjm,llm) ! pressure |
---|
| 135 | REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1) |
---|
| 136 | REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1) |
---|
| 137 | REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm) |
---|
| 138 | |
---|
| 139 | REAL(rstd) :: berni(iim*jjm,llm) ! (W/m_il)^2 |
---|
| 140 | REAL(rstd) :: berni1(iim*jjm) ! (W/m_il)^2 |
---|
| 141 | !REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff |
---|
| 142 | REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Cvd, Rd_preff |
---|
| 143 | INTEGER :: ij, l |
---|
| 144 | |
---|
| 145 | CALL trace_start("compute_caldyn_solver") |
---|
| 146 | |
---|
| 147 | !Rd=cpp*kappa |
---|
| 148 | |
---|
| 149 | IF(dysl) THEN |
---|
| 150 | |
---|
| 151 | !$OMP BARRIER |
---|
| 152 | |
---|
| 153 | #include "../kernels_hex/caldyn_mil.k90" |
---|
| 154 | IF(tau>0) THEN ! solve implicit problem for geopotential |
---|
| 155 | CALL compute_NH_geopot(tau,phis, rhodz, m_il, theta, W, geopot) |
---|
| 156 | END IF |
---|
| 157 | #define PHI_BOT(ij) phis(ij) |
---|
| 158 | #include "../kernels_hex/caldyn_solver.k90" |
---|
| 159 | #undef PHI_BOT |
---|
| 160 | !$OMP BARRIER |
---|
| 161 | |
---|
| 162 | ELSE |
---|
| 163 | |
---|
| 164 | #define BERNI(ij) berni1(ij) |
---|
| 165 | ! FIXME : vertical OpenMP parallelism will not work |
---|
| 166 | |
---|
| 167 | ! average m_ik to interfaces => m_il |
---|
| 168 | !DIR$ SIMD |
---|
| 169 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 170 | m_il(ij,1) = .5*rhodz(ij,1) |
---|
| 171 | ENDDO |
---|
| 172 | DO l=2,llm |
---|
| 173 | !DIR$ SIMD |
---|
| 174 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 175 | m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l)) |
---|
| 176 | ENDDO |
---|
| 177 | ENDDO |
---|
| 178 | !DIR$ SIMD |
---|
| 179 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 180 | m_il(ij,llm+1) = .5*rhodz(ij,llm) |
---|
| 181 | ENDDO |
---|
| 182 | |
---|
| 183 | IF(tau>0) THEN ! solve implicit problem for geopotential |
---|
| 184 | CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot) |
---|
| 185 | END IF |
---|
| 186 | |
---|
| 187 | ! Compute pressure, stored temporarily in pk |
---|
| 188 | ! kappa = R/Cp |
---|
| 189 | ! 1-kappa = Cv/Cp |
---|
| 190 | ! Cp/Cv = 1/(1-kappa) |
---|
| 191 | gamma = 1./(1.-kappa) |
---|
| 192 | DO l=1,llm |
---|
| 193 | !DIR$ SIMD |
---|
| 194 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 195 | rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l)) |
---|
| 196 | X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij |
---|
| 197 | ! kappa.theta.rho = p/exner |
---|
| 198 | ! => X = (p/p0)/(exner/Cp) |
---|
| 199 | ! = (p/p0)^(1-kappa) |
---|
| 200 | pk(ij,l) = preff*(X_ij**gamma) |
---|
| 201 | ENDDO |
---|
| 202 | ENDDO |
---|
| 203 | |
---|
| 204 | ! Update W, compute tendencies |
---|
| 205 | DO l=2,llm |
---|
| 206 | !DIR$ SIMD |
---|
| 207 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 208 | dW(ij,l) = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l) |
---|
| 209 | W(ij,l) = W(ij,l)+tau*dW(ij,l) ! update W |
---|
| 210 | dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l) |
---|
| 211 | ENDDO |
---|
| 212 | ! PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l))) |
---|
| 213 | ! PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l))) |
---|
| 214 | ENDDO |
---|
| 215 | ! Lower BC (FIXME : no orography yet !) |
---|
| 216 | DO ij=ij_begin,ij_end |
---|
| 217 | dPhi(ij,1)=0 |
---|
| 218 | W(ij,1)=0 |
---|
| 219 | dW(ij,1)=0 |
---|
| 220 | dPhi(ij,llm+1)=0 ! rigid lid |
---|
| 221 | W(ij,llm+1)=0 |
---|
| 222 | dW(ij,llm+1)=0 |
---|
| 223 | ENDDO |
---|
| 224 | ! Upper BC p=ptop |
---|
| 225 | ! DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 226 | ! dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm) |
---|
| 227 | ! dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm) |
---|
| 228 | ! ENDDO |
---|
| 229 | |
---|
| 230 | ! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2) |
---|
| 231 | DO l=1,llm |
---|
| 232 | !DIR$ SIMD |
---|
| 233 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 234 | pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow |
---|
| 235 | BERNI(ij) = (-.25*g*g)*( & |
---|
| 236 | (W(ij,l)/m_il(ij,l))**2 & |
---|
| 237 | + (W(ij,l+1)/m_il(ij,l+1))**2 ) |
---|
| 238 | ENDDO |
---|
| 239 | DO ij=ij_begin,ij_end |
---|
| 240 | du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right)) |
---|
| 241 | du(ij+u_lup,l) = ne_lup *(BERNI(ij)-BERNI(ij+t_lup)) |
---|
| 242 | du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) |
---|
| 243 | ENDDO |
---|
| 244 | ENDDO |
---|
| 245 | #undef BERNI |
---|
| 246 | |
---|
| 247 | END IF ! dysl |
---|
| 248 | |
---|
| 249 | CALL trace_end("compute_caldyn_solver") |
---|
| 250 | |
---|
| 251 | END SUBROUTINE compute_caldyn_solver |
---|
| 252 | |
---|
| 253 | END MODULE compute_caldyn_solver_mod |
---|