Changeset 878 for codes/icosagcm/devel/Python/src
- Timestamp:
- 05/29/19 20:33:00 (5 years ago)
- Location:
- codes/icosagcm/devel/Python/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin
r876 r878 10 10 {%- endmacro %} 11 11 12 KERNEL(compute_NH_geopot)13 tau2_g=tau*tau/g14 g2=g*g15 gm2 = 1./g216 vreff = Treff*cpp/preff*kappa17 gamma = 1./(1.-kappa)18 19 BARRIER20 21 ! compute Phi_star22 SEQUENCE_C123 BODY('1,llm+1')24 Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau)25 END_BLOCK26 END_BLOCK27 28 ! Newton-Raphson iteration : Phi_il contains current guess value29 DO iter=1,2 ! 2 iterations should be enough30 ! Compute pressure, A_ik31 SELECT CASE(caldyn_thermo)32 CASE(thermo_theta)33 {% call() compute_p_and_Aik() %}34 X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij35 p_ik(CELL) = preff*(X_ij**gamma)36 c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho37 {% endcall %}38 CASE(thermo_entropy)39 {% call() compute_p_and_Aik() %}40 X_ij = log(vreff*rho_ij) + theta(CELL)/cpp41 p_ik(CELL) = preff*exp(X_ij*gamma)42 c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho43 {% endcall %}44 CASE DEFAULT45 PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo46 STOP47 END SELECT48 49 ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom50 ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm51 52 SEQUENCE_C153 ! Compute residual R_il and B_il54 PROLOGUE(1)55 ! bottom interface l=156 ml_g2 = gm2*m_il(CELL)57 B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot58 R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL)) &59 + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) )60 END_BLOCK61 BODY('2,llm')62 ! inner interfaces63 ml_g2 = gm2*m_il(CELL)64 B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g265 R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL)) &66 + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL)))67 ! consistency check : if Wil=0 and initial state is in hydrostatic balance68 ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^269 ! and residual = tau^2*(ml+(1/g)dl_pi)=070 END_BLOCK71 EPILOGUE('llm+1')72 ! top interface l=llm+173 ml_g2 = gm2*m_il(CELL)74 B_il(CELL) = A_ik(DOWN(CELL)) + ml_g275 R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL)) &76 + tau2_g*( ptop-p_ik(DOWN(CELL)) )77 END_BLOCK78 !79 ! Forward sweep :80 ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),81 ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))82 PROLOGUE(1)83 X_ij = 1./B_il(CELL)84 C_ik(CELL) = -A_ik(CELL) * X_ij85 D_il(CELL) = R_il(CELL) * X_ij86 END_BLOCK87 BODY('2,llm')88 X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )89 C_ik(CELL) = -A_ik(CELL) * X_ij90 D_il(CELL) = (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij91 END_BLOCK92 EPILOGUE('llm+1')93 X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )94 D_il(CELL) = (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij95 ! Back substitution :96 ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=097 ! + Newton-Raphson update98 ! top interface l=llm+199 x_il(CELL) = D_il(CELL)100 Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)101 END_BLOCK102 BODY('llm,1,-1')103 ! Back substitution at lower interfaces104 x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL))105 Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)106 END_BLOCK107 END_BLOCK108 109 IF(debug_hevi_solver) THEN110 PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))111 PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))112 DO l=1,llm+1113 WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:)))114 END DO115 DO l=1,llm+1116 WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:)))117 END DO118 END IF119 120 END DO ! Newton-Raphson121 122 BARRIER123 124 debug_hevi_solver=.FALSE.125 END_BLOCK126 127 12 KERNEL(caldyn_mil) 128 13 FORALL_CELLS_EXT('1', 'llm+1') … … 131 16 CST_IF(IS_INNER_INTERFACE, m_il(CELL) = .5*(rhodz(CELL)+rhodz(DOWN(CELL))) ) 132 17 CST_IF(IS_TOP_INTERFACE, m_il(CELL) = .5*rhodz(DOWN(CELL)) ) 133 END_BLOCK134 END_BLOCK135 END_BLOCK136 137 KERNEL(caldyn_solver)138 !139 ! Compute pressure (pres) and Exner function (pk)140 ! kappa = R/Cp141 ! 1-kappa = Cv/Cp142 ! Cp/Cv = 1/(1-kappa)143 gamma = 1./(1.-kappa)144 vreff = Rd*Treff/preff ! reference specific volume145 Cvd = 1./(cpp-Rd)146 Rd_preff = kappa*cpp/preff147 FORALL_CELLS_EXT()148 SELECT CASE(caldyn_thermo)149 CASE(thermo_theta)150 ON_PRIMAL151 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_ij154 ! kappa.theta.rho = p/exner155 ! => X = (p/p0)/(exner/Cp)156 ! = (p/p0)^(1-kappa)157 pres(CELL) = preff*(X_ij**gamma) ! pressure158 ! Compute Exner function (needed by compute_caldyn_fast)159 ! other formulae possible if exponentiation is slow160 pk(CELL) = cpp*((pres(CELL)/preff)**kappa) ! Exner161 END_BLOCK162 CASE(thermo_entropy)163 ON_PRIMAL164 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 )167 pres(CELL) = rho_ij*Rd*T_ij168 pk(CELL) = T_ij169 END_BLOCK170 CASE DEFAULT171 STOP172 END SELECT173 END_BLOCK174 175 ! We need a barrier here because we compute pres above and do a vertical difference below176 BARRIER177 178 FORALL_CELLS_EXT('1', 'llm+1')179 ON_PRIMAL180 CST_IFTHEN(IS_BOTTOM_LEVEL)181 ! Lower BC182 dW(CELL) = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL)183 CST_ELSEIF(IS_TOP_INTERFACE)184 ! Top BC185 dW(CELL) = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL)186 CST_ELSE187 dW(CELL) = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL)188 CST_ENDIF189 W(CELL) = W(CELL)+tau*dW(CELL) ! update W190 dPhi(CELL) = g*g*W(CELL)/m_il(CELL)191 END_BLOCK192 END_BLOCK193 194 ! We need a barrier here because we update W above and do a vertical average below195 BARRIER196 197 FORALL_CELLS_EXT()198 ON_PRIMAL199 ! compute du = -0.5*g^2.grad(w^2)200 berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 )201 END_BLOCK202 END_BLOCK203 FORALL_CELLS_EXT()204 ON_EDGES205 du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))206 18 END_BLOCK 207 19 END_BLOCK -
codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin
r876 r878 96 96 END_BLOCK 97 97 98 KERNEL(theta)99 IF(caldyn_eta==eta_mass) THEN ! Compute mass100 ! FIXME : here mass_col is computed from rhodz101 ! so that the DOFs are the same whatever caldyn_eta102 ! in DYNAMICO mass_col is prognosed rather than rhodz103 SEQUENCE_C1104 PROLOGUE(0)105 mass_col(HIDX(CELL))=0.106 END_BLOCK107 BODY('1,llm')108 mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL)109 END_BLOCK110 END_BLOCK111 FORALL_CELLS_EXT()112 ON_PRIMAL113 ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass114 ! m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL)115 ! rhodz(CELL) = m/g116 rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL)117 END_BLOCK118 END_BLOCK119 END IF120 DO iq=1,nqdyn121 FORALL_CELLS_EXT()122 ON_PRIMAL123 theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL)124 END_BLOCK125 END_BLOCK126 END DO127 END_BLOCK128 98
Note: See TracChangeset
for help on using the changeset viewer.