Changeset 878 for codes/icosagcm/devel/src
- Timestamp:
- 05/29/19 20:33:00 (5 years ago)
- Location:
- codes/icosagcm/devel/src
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dynamics/compute_NH_geopot.F90
r859 r878 6 6 LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 7 7 8 PUBLIC :: compute_NH_geopot 8 #include "../unstructured/unstructured.h90" 9 10 PUBLIC :: compute_NH_geopot,compute_NH_geopot_unst 9 11 10 12 CONTAINS 13 14 #ifdef BEGIN_DYSL 15 16 KERNEL(compute_NH_geopot) 17 tau2_g=tau*tau/g 18 g2=g*g 19 gm2 = 1./g2 20 vreff = Treff*cpp/preff*kappa 21 gamma = 1./(1.-kappa) 22 23 BARRIER 24 25 ! compute Phi_star 26 SEQUENCE_C1 27 BODY('1,llm+1') 28 Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau) 29 END_BLOCK 30 END_BLOCK 31 32 ! Newton-Raphson iteration : Phi_il contains current guess value 33 DO iter=1,2 ! 2 iterations should be enough 34 ! Compute pressure, A_ik 35 SELECT CASE(caldyn_thermo) 36 CASE(thermo_theta) 37 {% call() compute_p_and_Aik() %} 38 X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij 39 p_ik(CELL) = preff*(X_ij**gamma) 40 c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho 41 {% endcall %} 42 CASE(thermo_entropy) 43 {% call() compute_p_and_Aik() %} 44 X_ij = log(vreff*rho_ij) + theta(CELL)/cpp 45 p_ik(CELL) = preff*exp(X_ij*gamma) 46 c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho 47 {% endcall %} 48 CASE DEFAULT 49 PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo 50 STOP 51 END SELECT 52 53 ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom 54 ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm 55 56 SEQUENCE_C1 57 ! Compute residual R_il and B_il 58 PROLOGUE(1) 59 ! bottom interface l=1 60 ml_g2 = gm2*m_il(CELL) 61 B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot 62 R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL)) & 63 + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) ) 64 END_BLOCK 65 BODY('2,llm') 66 ! inner interfaces 67 ml_g2 = gm2*m_il(CELL) 68 B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g2 69 R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL)) & 70 + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL))) 71 ! consistency check : if Wil=0 and initial state is in hydrostatic balance 72 ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^2 73 ! and residual = tau^2*(ml+(1/g)dl_pi)=0 74 END_BLOCK 75 EPILOGUE('llm+1') 76 ! top interface l=llm+1 77 ml_g2 = gm2*m_il(CELL) 78 B_il(CELL) = A_ik(DOWN(CELL)) + ml_g2 79 R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL)) & 80 + tau2_g*( ptop-p_ik(DOWN(CELL)) ) 81 END_BLOCK 82 ! 83 ! Forward sweep : 84 ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)), 85 ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1)) 86 PROLOGUE(1) 87 X_ij = 1./B_il(CELL) 88 C_ik(CELL) = -A_ik(CELL) * X_ij 89 D_il(CELL) = R_il(CELL) * X_ij 90 END_BLOCK 91 BODY('2,llm') 92 X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) ) 93 C_ik(CELL) = -A_ik(CELL) * X_ij 94 D_il(CELL) = (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij 95 END_BLOCK 96 EPILOGUE('llm+1') 97 X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) ) 98 D_il(CELL) = (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij 99 ! Back substitution : 100 ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0 101 ! + Newton-Raphson update 102 ! top interface l=llm+1 103 x_il(CELL) = D_il(CELL) 104 Phi_il(CELL) = Phi_il(CELL) - x_il(CELL) 105 END_BLOCK 106 BODY('llm,1,-1') 107 ! Back substitution at lower interfaces 108 x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL)) 109 Phi_il(CELL) = Phi_il(CELL) - x_il(CELL) 110 END_BLOCK 111 END_BLOCK 112 113 IF(debug_hevi_solver) THEN 114 PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il)) 115 PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il)) 116 DO l=1,llm+1 117 WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:))) 118 END DO 119 DO l=1,llm+1 120 WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:))) 121 END DO 122 END IF 123 124 END DO ! Newton-Raphson 125 126 BARRIER 127 128 debug_hevi_solver=.FALSE. 129 END_BLOCK 130 131 #endif END_DYSL 132 133 SUBROUTINE compute_NH_geopot_unst(tau, m_ik, m_il, theta, W_il, Phi_il) 134 USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT 135 USE disvert_mod, only : g,Treff,cpp,preff,kappa,caldyn_thermo,thermo_theta, & 136 thermo_entropy 137 USE disvert_mod, ONLY : ptop 138 USE data_unstructured_mod, ONLY : primal_num,edge_num,dual_num,id_NH_geopot,debug_hevi_solver_, & 139 PHI_BOT,pbot,rho_bot,enter_trace, exit_trace 140 FIELD_MASS :: m_ik, theta ! IN*2 141 FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il ! IN,INOUT*2, LOCAL 142 NUM :: tau, gamma, tau2_g, tau2_g2, g2, gm2, vreff, Rd_preff 143 INTEGER :: iter 144 LOGICAL :: debug_hevi_solver 145 DECLARE_INDICES 146 NUM :: rho_ij, X_ij, Y_ij, wil, rho_c2_mik, c2_mik, ml_g2 147 #define COLUMN 0 148 #if COLUMN 149 NUM1(llm) :: pk, Ak, Ck 150 NUM1(llm+1):: Rl, Bl, Dl, xl 151 #define p_ik(l,ij) pk(l) 152 #define A_ik(l,ij) Ak(l) 153 #define C_ik(l,ij) Ck(l) 154 #define R_il(l,ij) Rl(l) 155 #define B_il(l,ij) Bl(l) 156 #define D_il(l,ij) Dl(l) 157 #define x_il(l,ij) xl(l) 158 #else 159 FIELD_MASS :: p_ik, A_ik, C_ik 160 FIELD_GEOPOT :: R_il, B_il, D_il, x_il 161 #endif 162 163 debug_hevi_solver=.FALSE. 164 !$OMP MASTER 165 debug_hevi_solver = debug_hevi_solver_ 166 !$OMP END MASTER 167 168 #define PHI_BOT(ij) Phi_bot 169 START_TRACE(id_NH_geopot, 7,0,0) 170 #include "../kernels_unst/compute_NH_geopot.k90" 171 STOP_TRACE 172 173 !$OMP MASTER 174 debug_hevi_solver_ = debug_hevi_solver 175 !$OMP END MASTER 176 #undef PHI_BOT 177 178 #if COLUMN 179 #undef p_ik 180 #undef A_ik 181 #undef C_ik 182 #undef R_il 183 #undef B_il 184 #undef D_il 185 #undef x_il 186 #endif 187 #undef COLUMN 188 END SUBROUTINE compute_NH_geopot_unst 11 189 12 190 SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il) -
codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90
r859 r878 4 4 PRIVATE 5 5 6 #include "../unstructured/unstructured.h90" 7 6 8 PUBLIC :: compute_caldyn_solver 7 9 8 10 CONTAINS 11 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 9 116 10 117 SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du) -
codes/icosagcm/devel/src/dynamics/compute_theta.F90
r831 r878 4 4 PRIVATE 5 5 6 #include "../unstructured/unstructured.h90" 7 6 8 PUBLIC :: compute_theta 7 9 8 10 CONTAINS 11 12 #ifdef BEGIN_DYSL 13 14 KERNEL(theta) 15 IF(caldyn_eta==eta_mass) THEN ! Compute mass 16 ! FIXME : here mass_col is computed from rhodz 17 ! so that the DOFs are the same whatever caldyn_eta 18 ! in DYNAMICO mass_col is prognosed rather than rhodz 19 SEQUENCE_C1 20 PROLOGUE(0) 21 mass_col(HIDX(CELL))=0. 22 END_BLOCK 23 BODY('1,llm') 24 mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL) 25 END_BLOCK 26 END_BLOCK 27 FORALL_CELLS_EXT() 28 ON_PRIMAL 29 ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on 30 ! pressure rather than mass 31 ! m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL) 32 ! rhodz(CELL) = m/g 33 rhodz(CELL) = MASS_DAK(CELL) + mass_col(HIDX(CELL))*MASS_DBK(CELL) 34 END_BLOCK 35 END_BLOCK 36 END IF 37 DO iq=1,nqdyn 38 FORALL_CELLS_EXT() 39 ON_PRIMAL 40 theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL) 41 END_BLOCK 42 END_BLOCK 43 END DO 44 END_BLOCK 45 46 #endif END_DYSL 47 48 SUBROUTINE compute_theta_unst(mass_col,rhodz,theta_rhodz, theta) 49 USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT 50 USE data_unstructured_mod, ONLY : id_theta,primal_num,dual_num,edge_num, & 51 enter_trace, exit_trace 52 USE grid_param, ONLY : nqdyn 53 USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop 54 FIELD_PS :: mass_col 55 FIELD_MASS :: rhodz 56 FIELD_THETA :: theta, theta_rhodz 57 DECLARE_INDICES 58 NUM :: m 59 START_TRACE(id_theta, 3,0,0) ! primal, dual, edge 60 #define MASS_DAK(l,ij) mass_dak(l) 61 #define MASS_DBK(l,ij) mass_dbk(l) 62 #include "../kernels_unst/theta.k90" 63 #undef MASS_DAK 64 #undef MASS_DBK 65 STOP_TRACE 66 END SUBROUTINE 9 67 10 68 SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta) -
codes/icosagcm/devel/src/kernels_unst/theta.k90
r686 r878 20 20 ! m = mass_dak(l,ij)+(mass_col(ij)*g+ptop)*mass_dbk(l,ij) 21 21 ! rhodz(l,ij) = m/g 22 rhodz(l,ij) = mass_dak(l,ij) + mass_col(ij)*mass_dbk(l,ij)22 rhodz(l,ij) = MASS_DAK(l,ij) + mass_col(ij)*MASS_DBK(l,ij) 23 23 END DO 24 24 END DO
Note: See TracChangeset
for help on using the changeset viewer.