Changeset 659 for codes/icosagcm/devel
- Timestamp:
- 12/30/17 02:01:51 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin
r615 r659 14 14 g2=g*g 15 15 gm2 = 1./g2 16 vreff = Treff*cpp/preff*kappa 16 17 gamma = 1./(1.-kappa) 17 18 … … 37 38 CASE(thermo_entropy) 38 39 {% call() compute_p_and_Aik() %} 39 X_ij = Treff*exp(theta(CELL)/cpp) ! theta = Tref.exp(s/Cp) 40 X_ij = (cpp/preff)*kappa*X_ij*rho_ij 41 p_ik(CELL) = preff*(X_ij**gamma) 40 X_ij = log(vreff*rho_ij) + theta(CELL)/cpp 41 p_ik(CELL) = preff*exp(X_ij*gamma) 42 42 c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho 43 43 {% endcall %} … … 125 125 END_BLOCK 126 126 127 KERNEL(caldyn_ solver)127 KERNEL(caldyn_mil) 128 128 FORALL_CELLS_EXT('1', 'llm+1') 129 129 ON_PRIMAL … … 133 133 END_BLOCK 134 134 END_BLOCK 135 ! 136 IF(tau>0) THEN ! solve implicit problem for geopotential 137 CALL compute_NH_geopot(tau,PHI_BOT_VAR, rhodz, m_il, theta, W, geopot) 138 END IF 135 END_BLOCK 136 137 KERNEL(caldyn_solver) 139 138 ! 140 139 ! Compute pressure (pres) and Exner function (pk) … … 142 141 ! 1-kappa = Cv/Cp 143 142 ! Cp/Cv = 1/(1-kappa) 144 gamma = 1./(1.-kappa)143 gamma = 1./(1.-kappa) 145 144 vreff = Rd*Treff/preff ! reference specific volume 146 Cvd = cpp-Rd 145 Cvd = 1./(cpp-Rd) 146 Rd_preff = kappa*cpp/preff 147 147 FORALL_CELLS_EXT() 148 148 SELECT CASE(caldyn_thermo) 149 149 CASE(thermo_theta) 150 150 ON_PRIMAL 151 rho_ij = g*rhodz(CELL)/(geopot(UP(CELL))-geopot(CELL)) 152 X_ij = (cpp/preff)*kappa*theta(CELL,1)*rho_ij 151 rho_ij = 1./(geopot(UP(CELL))-geopot(CELL)) 152 rho_ij = rho_ij*g*rhodz(CELL) 153 X_ij = Rd_preff*theta(CELL,1)*rho_ij 153 154 ! kappa.theta.rho = p/exner 154 155 ! => X = (p/p0)/(exner/Cp) … … 161 162 CASE(thermo_entropy) 162 163 ON_PRIMAL 163 rho_ij = g*rhodz(CELL)/(geopot(UP(CELL))-geopot(CELL)) 164 T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))/Cvd ) 164 rho_ij = 1./(geopot(UP(CELL))-geopot(CELL)) 165 rho_ij = rho_ij*g*rhodz(CELL) 166 T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd ) 165 167 pres(CELL) = rho_ij*Rd*T_ij 166 168 pk(CELL) = T_ij
Note: See TracChangeset
for help on using the changeset viewer.