[920] | 1 | MODULE compute_geopot_mod |
---|
| 2 | USE prec, ONLY : rstd |
---|
| 3 | USE grid_param |
---|
| 4 | USE earth_const |
---|
| 5 | USE disvert_mod |
---|
| 6 | USE omp_para |
---|
| 7 | USE trace |
---|
| 8 | IMPLICIT NONE |
---|
| 9 | PRIVATE |
---|
| 10 | SAVE |
---|
| 11 | |
---|
[937] | 12 | #include "../unstructured/unstructured.h90" |
---|
[920] | 13 | |
---|
[937] | 14 | PUBLIC :: compute_geopot_unst, compute_geopot_hex, compute_geopot_manual |
---|
| 15 | |
---|
[920] | 16 | CONTAINS |
---|
| 17 | |
---|
[937] | 18 | #ifdef BEGIN_DYSL |
---|
| 19 | |
---|
| 20 | #define THECELL {{ thecell }} |
---|
| 21 | |
---|
| 22 | {# ---------------- macro to generate code computing pressure top-down --------------- |
---|
| 23 | formula = formula to compute 'gravitational' mass |
---|
| 24 | = rhodz (dry) rhodz*theta (boussinesq) rhodz*(1+qv) (moist) #} |
---|
| 25 | |
---|
| 26 | #define BALANCE(formula) {% call(thecell) balance() %} formula {% endcall %} |
---|
| 27 | {% macro balance() %} {% set formula=caller %} |
---|
| 28 | SEQUENCE_C1 |
---|
| 29 | PROLOGUE(llm) |
---|
| 30 | pk(CELL) = ptop + .5*g*{{ formula('CELL') }} |
---|
| 31 | END_BLOCK |
---|
| 32 | BODY('llm-1,1,-1') |
---|
| 33 | pk(CELL) = pk(UP(CELL)) + (.5*g)*({{ formula('CELL') }}+{{ formula('UP(CELL)') }}) |
---|
| 34 | END_BLOCK |
---|
| 35 | IF(caldyn_eta == eta_lag) THEN |
---|
| 36 | EPILOGUE(1) |
---|
| 37 | ps(HIDX(CELL)) = pk(CELL) + .5*g*{{ formula('CELL') }} |
---|
| 38 | END_BLOCK |
---|
| 39 | END IF |
---|
| 40 | END_BLOCK |
---|
| 41 | {%- endmacro %} |
---|
| 42 | |
---|
| 43 | {# ------------ macro to generate code computing geopotential bottom-up -------------- |
---|
| 44 | var = variable to be stored in pk(CELL) |
---|
| 45 | caller() computes gv = g*v where v = specific volume |
---|
| 46 | details depend on caldyn_thermo #} |
---|
| 47 | |
---|
| 48 | #define GEOPOT(var) {% call geopot(var) %} |
---|
| 49 | {% macro geopot(var) %} {% set formula=caller %} |
---|
| 50 | SEQUENCE_C1 |
---|
| 51 | BODY('1,llm') |
---|
| 52 | p_ik = pk(CELL) |
---|
| 53 | {{ formula() }} |
---|
| 54 | pk(CELL) = {{ var }} |
---|
| 55 | geopot(UP(CELL)) = geopot(CELL) + gv*rhodz(CELL) |
---|
| 56 | END_BLOCK |
---|
| 57 | END_BLOCK |
---|
| 58 | {%- endmacro %} |
---|
| 59 | |
---|
| 60 | #define END_GEOPOT {% endcall %} |
---|
| 61 | |
---|
| 62 | KERNEL(compute_geopot) |
---|
| 63 | SELECT CASE(caldyn_thermo) |
---|
| 64 | CASE(thermo_boussinesq) |
---|
| 65 | ! use hydrostatic balance with theta*rhodz to find pk (=Lagrange multiplier=pressure) |
---|
| 66 | BALANCE( theta(THECELL,1)*rhodz(THECELL) ) |
---|
| 67 | ! now pk contains the Lagrange multiplier (pressure) |
---|
| 68 | ! specific volume 1 = dphi/g/rhodz |
---|
| 69 | SEQUENCE_C1 |
---|
| 70 | BODY('1,llm') |
---|
| 71 | geopot(UP(CELL)) = geopot(CELL) + g*rhodz(CELL) |
---|
| 72 | END_BLOCK |
---|
| 73 | END_BLOCK |
---|
| 74 | CASE(thermo_theta) |
---|
| 75 | BALANCE( rhodz(THECELL) ) |
---|
| 76 | GEOPOT('exner_ik') |
---|
| 77 | exner_ik = cpp * (p_ik/preff) ** kappa |
---|
| 78 | gv = (g*kappa)*theta(CELL,1)*exner_ik/p_ik |
---|
| 79 | END_GEOPOT |
---|
| 80 | CASE(thermo_entropy) |
---|
| 81 | BALANCE( rhodz(THECELL) ) |
---|
| 82 | GEOPOT('temp_ik') |
---|
| 83 | temp_ik = Treff*exp((theta(CELL,1) + Rd*log(p_ik/preff))/cpp) |
---|
| 84 | gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p |
---|
| 85 | END_GEOPOT |
---|
| 86 | CASE(thermo_variable_Cp) |
---|
| 87 | ! thermodynamics with variable Cp |
---|
| 88 | ! Cp.dT = dh = Tds + vdp |
---|
| 89 | ! pv = RT |
---|
| 90 | ! => ds = (dh+v.dp)/T = Cp.dT/T - R dp/p |
---|
| 91 | ! Cp(T) = Cp0 * (T/T0)^nu |
---|
| 92 | ! => s(p,T) = Cp(T)/nu - R log(p/preff) |
---|
| 93 | ! h = Cp(T).T/(nu+1) |
---|
| 94 | BALANCE( rhodz(THECELL) ) |
---|
| 95 | GEOPOT('temp_ik') |
---|
| 96 | Cp_ik = nu*( theta(CELL,1) + Rd*log(p_ik/preff) ) |
---|
| 97 | temp_ik = Treff* (Cp_ik/cpp)**(1./nu) |
---|
| 98 | gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p |
---|
| 99 | END_GEOPOT |
---|
| 100 | CASE(thermo_moist) |
---|
| 101 | BALANCE( rhodz(THECELL)*(1.+theta(THECELL,2)) ) |
---|
| 102 | GEOPOT('temp_ik') |
---|
| 103 | qv = theta(CELL,2) ! water vaper mixing ratio = mv/md |
---|
| 104 | Rmix = Rd+qv*Rv |
---|
| 105 | chi = ( theta(CELL,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff) |
---|
| 106 | temp_ik = Treff*exp(chi) |
---|
| 107 | ! specific volume v = R*T/p |
---|
| 108 | ! R = (Rd + qv.Rv)/(1+qv) |
---|
| 109 | gv = g*Rmix*temp_ik/(p_ik*(1+qv)) |
---|
| 110 | END_GEOPOT |
---|
| 111 | END SELECT |
---|
| 112 | END_BLOCK |
---|
| 113 | |
---|
| 114 | #endif END_DYSL |
---|
| 115 | |
---|
[920] | 116 | !**************************** Geopotential ***************************** |
---|
| 117 | |
---|
[937] | 118 | SUBROUTINE compute_geopot_unst(rhodz,theta,ps,pk,geopot) |
---|
| 119 | USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT |
---|
| 120 | USE data_unstructured_mod, ONLY : enter_trace, exit_trace, & |
---|
| 121 | id_geopot |
---|
| 122 | |
---|
| 123 | FIELD_MASS :: rhodz,pk ! IN, OUT |
---|
| 124 | FIELD_THETA :: theta ! IN |
---|
| 125 | FIELD_GEOPOT :: geopot ! IN(l=1)/OUT(l>1) |
---|
| 126 | FIELD_PS :: ps ! OUT |
---|
| 127 | DECLARE_INDICES |
---|
| 128 | NUM :: chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix, Cp_ik |
---|
| 129 | START_TRACE(id_geopot, 3,0,3) |
---|
| 130 | #include "../kernels_unst/compute_geopot.k90" |
---|
| 131 | STOP_TRACE |
---|
| 132 | END SUBROUTINE compute_geopot_unst |
---|
| 133 | |
---|
[920] | 134 | SUBROUTINE compute_geopot_hex(rhodz,theta, ps,pk,geopot) |
---|
| 135 | REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) |
---|
| 136 | REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn) ! active scalars : theta/entropy, moisture, ... |
---|
| 137 | REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) |
---|
| 138 | REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm) ! Exner function (compressible) /Lagrange multiplier (Boussinesq) |
---|
| 139 | REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential |
---|
| 140 | |
---|
| 141 | INTEGER :: i,j,ij,l |
---|
| 142 | REAL(rstd) :: p_ik, exner_ik, Cp_ik, temp_ik, qv, chi, Rmix, gv |
---|
| 143 | INTEGER :: ij_omp_begin_ext, ij_omp_end_ext |
---|
| 144 | |
---|
| 145 | CALL trace_start("compute_geopot") |
---|
| 146 | !$OMP BARRIER |
---|
| 147 | CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext) |
---|
| 148 | #include "../kernels_hex/compute_geopot.k90" |
---|
| 149 | !$OMP BARRIER |
---|
| 150 | |
---|
| 151 | CALL trace_end("compute_geopot") |
---|
| 152 | END SUBROUTINE compute_geopot_hex |
---|
| 153 | |
---|
| 154 | SUBROUTINE compute_geopot_manual(rhodz,theta, ps,pk,geopot) |
---|
| 155 | REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) |
---|
| 156 | REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn) ! active scalars : theta/entropy, moisture, ... |
---|
| 157 | REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) |
---|
| 158 | REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm) ! Exner function (compressible) /Lagrange multiplier (Boussinesq) |
---|
| 159 | REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential |
---|
| 160 | |
---|
| 161 | INTEGER :: i,j,ij,l |
---|
| 162 | REAL(rstd) :: p_ik, exner_ik, Cp_ik, temp_ik, qv, chi, Rmix, gv |
---|
| 163 | INTEGER :: ij_omp_begin_ext, ij_omp_end_ext |
---|
| 164 | |
---|
| 165 | CALL trace_start("compute_geopot") |
---|
| 166 | |
---|
| 167 | !$OMP BARRIER |
---|
| 168 | |
---|
| 169 | CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext) |
---|
| 170 | |
---|
| 171 | ! Pressure is computed first top-down (temporarily stored in pk) |
---|
| 172 | ! Then Exner pressure and geopotential are computed bottom-up |
---|
| 173 | ! Works also when caldyn_eta=eta_mass |
---|
| 174 | |
---|
| 175 | IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier |
---|
| 176 | ! specific volume 1 = dphi/g/rhodz |
---|
| 177 | ! IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency |
---|
| 178 | DO l = 1,llm |
---|
| 179 | !DIR$ SIMD |
---|
| 180 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 181 | geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) |
---|
| 182 | ENDDO |
---|
| 183 | ENDDO |
---|
| 184 | ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) |
---|
| 185 | ! uppermost layer |
---|
| 186 | !DIR$ SIMD |
---|
| 187 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 188 | pk(ij,llm) = ptop + (.5*g)*theta(ij,llm,1)*rhodz(ij,llm) |
---|
| 189 | END DO |
---|
| 190 | ! other layers |
---|
| 191 | DO l = llm-1, 1, -1 |
---|
| 192 | ! !$OMP DO SCHEDULE(STATIC) |
---|
| 193 | !DIR$ SIMD |
---|
| 194 | DO ij=ij_begin_ext,ij_end_ext |
---|
| 195 | pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l,1)*rhodz(ij,l)+theta(ij,l+1,1)*rhodz(ij,l+1)) |
---|
| 196 | END DO |
---|
| 197 | END DO |
---|
| 198 | ! now pk contains the Lagrange multiplier (pressure) |
---|
| 199 | ELSE ! non-Boussinesq, compute pressure, Exner pressure or temperature, then geopotential |
---|
| 200 | ! uppermost layer |
---|
| 201 | |
---|
| 202 | SELECT CASE(caldyn_thermo) |
---|
| 203 | CASE(thermo_theta, thermo_entropy) |
---|
| 204 | !DIR$ SIMD |
---|
| 205 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 206 | pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) |
---|
| 207 | END DO |
---|
| 208 | ! other layers |
---|
| 209 | DO l = llm-1, 1, -1 |
---|
| 210 | !DIR$ SIMD |
---|
| 211 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 212 | pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) |
---|
| 213 | END DO |
---|
| 214 | END DO |
---|
| 215 | ! surface pressure (for diagnostics) |
---|
| 216 | IF(caldyn_eta==eta_lag) THEN |
---|
| 217 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 218 | ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) |
---|
| 219 | END DO |
---|
| 220 | END IF |
---|
| 221 | CASE(thermo_moist) ! theta(ij,l,2) = qv = mv/md |
---|
| 222 | !DIR$ SIMD |
---|
| 223 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 224 | pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)*(1.+theta(ij,l,2)) |
---|
| 225 | END DO |
---|
| 226 | ! other layers |
---|
| 227 | DO l = llm-1, 1, -1 |
---|
| 228 | !DIR$ SIMD |
---|
| 229 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 230 | pk(ij,l) = pk(ij,l+1) + (.5*g)*( & |
---|
| 231 | rhodz(ij,l) *(1.+theta(ij,l,2)) + & |
---|
| 232 | rhodz(ij,l+1)*(1.+theta(ij,l+1,2)) ) |
---|
| 233 | END DO |
---|
| 234 | END DO |
---|
| 235 | ! surface pressure (for diagnostics) |
---|
| 236 | IF(caldyn_eta==eta_lag) THEN |
---|
| 237 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 238 | ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)*(1.+theta(ij,l,2)) |
---|
| 239 | END DO |
---|
| 240 | END IF |
---|
| 241 | END SELECT |
---|
| 242 | |
---|
| 243 | DO l = 1,llm |
---|
| 244 | SELECT CASE(caldyn_thermo) |
---|
| 245 | CASE(thermo_theta) |
---|
| 246 | !DIR$ SIMD |
---|
| 247 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 248 | p_ik = pk(ij,l) |
---|
| 249 | exner_ik = cpp * (p_ik/preff) ** kappa |
---|
| 250 | pk(ij,l) = exner_ik |
---|
| 251 | ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz |
---|
| 252 | geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik |
---|
| 253 | ENDDO |
---|
| 254 | CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff) |
---|
| 255 | !DIR$ SIMD |
---|
| 256 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 257 | p_ik = pk(ij,l) |
---|
| 258 | temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp) |
---|
| 259 | pk(ij,l) = temp_ik |
---|
| 260 | ! specific volume v = Rd*T/p = dphi/g/rhodz |
---|
| 261 | geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik |
---|
| 262 | ENDDO |
---|
| 263 | CASE(thermo_moist) ! theta is moist pseudo-entropy per dry air mass |
---|
| 264 | DO ij=ij_omp_begin_ext,ij_omp_end_ext |
---|
| 265 | p_ik = pk(ij,l) |
---|
| 266 | qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md |
---|
| 267 | Rmix = Rd+qv*Rv |
---|
| 268 | chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff) |
---|
| 269 | temp_ik = Treff*exp(chi) |
---|
| 270 | pk(ij,l) = temp_ik |
---|
| 271 | ! specific volume v = R*T/p = dphi/g/rhodz |
---|
| 272 | ! R = (Rd + qv.Rv)/(1+qv) |
---|
| 273 | geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv)) |
---|
| 274 | ENDDO |
---|
| 275 | CASE DEFAULT |
---|
| 276 | STOP |
---|
| 277 | END SELECT |
---|
| 278 | ENDDO |
---|
| 279 | END IF |
---|
| 280 | |
---|
| 281 | !ym flush geopot |
---|
| 282 | !$OMP BARRIER |
---|
| 283 | |
---|
| 284 | CALL trace_end("compute_geopot") |
---|
| 285 | |
---|
| 286 | END SUBROUTINE compute_geopot_manual |
---|
| 287 | |
---|
| 288 | END MODULE compute_geopot_mod |
---|