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