Changeset 914 for codes/icosagcm/devel/src
- Timestamp:
- 06/17/19 18:05:22 (5 years ago)
- Location:
- codes/icosagcm/devel/src/diagnostics
- Files:
-
- 2 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/diagnostics/observable.f90
r913 r914 44 44 USE earth_const 45 45 USE compute_pression_mod, ONLY : pression_mid 46 USE compute_temperature_mod 47 USE compute_velocity_mod 46 48 USE vertical_interp_mod 47 49 USE theta2theta_rhodz_mod … … 97 99 98 100 CALL pression_mid(f_ps, f_pmid) 99 CALL diagnose_temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T101 CALL temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T 100 102 101 103 IF(init) THEN … … 125 127 IF(grid_type == grid_unst) RETURN 126 128 127 CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)129 CALL velocity(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_Fel, f_buf_uh, f_buf_i) 128 130 CALL transfert_request(f_buf_uh,req_e1_vect) 129 131 CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat) … … 210 212 END SUBROUTINE output_energyflux 211 213 212 !------------------- Conversion from prognostic to observable variables ------------------213 214 SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)215 USE disvert_mod, ONLY : caldyn_eta, eta_mass216 USE compute_diagnostics_mod, ONLY : compute_rhodz217 TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &218 f_u(:), f_W(:), f_uz(:), & ! IN219 f_uh(:) ! OUT220 REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:), F_el(:,:)221 INTEGER :: ind222 223 DO ind=1,ndomain224 IF (.NOT. assigned_domain(ind)) CYCLE225 CALL swap_dimensions(ind)226 CALL swap_geometry(ind)227 geopot = f_geopot(ind)228 rhodz = f_rhodz(ind)229 u = f_u(ind)230 W = f_W(ind)231 uh = f_uh(ind)232 F_el = f_buf_Fel(ind)233 IF(caldyn_eta==eta_mass) THEN234 ps=f_ps(ind)235 CALL compute_rhodz(.TRUE., ps, rhodz)236 END IF237 uz = f_uz(ind)238 !$OMP BARRIER239 CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W, F_el, uh,uz)240 !$OMP BARRIER241 END DO242 END SUBROUTINE progonostic_vel_to_horiz243 244 SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, F_el, uh, uz)245 USE omp_para246 REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)247 REAL(rstd), INTENT(IN) :: rhodz(iim*jjm,llm)248 REAL(rstd), INTENT(IN) :: u(3*iim*jjm,llm)249 REAL(rstd), INTENT(IN) :: W(iim*jjm,llm+1)250 REAL(rstd), INTENT(OUT) :: uh(3*iim*jjm,llm)251 REAL(rstd), INTENT(OUT) :: uz(iim*jjm,llm)252 INTEGER :: ij,l253 REAL(rstd) :: F_el(3*iim*jjm,llm+1)254 REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil255 ! NB : u and uh are not in DEC form, they are normal components256 ! => we must divide by de257 IF(hydrostatic) THEN258 uh(:,:)=u(:,:)259 uz(:,:)=0.260 ELSE261 DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)262 DO ij=ij_begin_ext, ij_end_ext263 ! Compute on edge 'right'264 W_el = .5*( W(ij,l)+W(ij+t_right,l) )265 DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))266 F_el(ij+u_right,l) = DePhil*W_el/de(ij+u_right)267 ! Compute on edge 'lup'268 W_el = .5*( W(ij,l)+W(ij+t_lup,l) )269 DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))270 F_el(ij+u_lup,l) = DePhil*W_el/de(ij+u_lup)271 ! Compute on edge 'ldown'272 W_el = .5*( W(ij,l)+W(ij+t_ldown,l) )273 DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))274 F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown)275 END DO276 END DO277 278 ! We need a barrier here because we compute F_el above and do a vertical average below279 !$OMP BARRIER280 281 DO l=ll_begin, ll_end ! compute on k levels (full levels)282 DO ij=ij_begin_ext, ij_end_ext283 ! w = vertical momentum = g^-2*dPhi/dt = uz/g284 uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)285 ! uh = u-w.grad(Phi) = u - uz.grad(z)286 uh(ij+u_right,l) = u(ij+u_right,l) - (F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) / (rhodz(ij,l)+rhodz(ij+t_right,l))287 uh(ij+u_lup,l) = u(ij+u_lup,l) - (F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1)) / (rhodz(ij,l)+rhodz(ij+t_lup,l))288 uh(ij+u_ldown,l) = u(ij+u_ldown,l) - (F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))289 END DO290 END DO291 292 END IF293 END SUBROUTINE compute_prognostic_vel_to_horiz294 295 SUBROUTINE diagnose_temperature(f_pmid,f_q,f_temp)296 USE icosa297 IMPLICIT NONE298 TYPE(t_field), POINTER :: f_pmid(:) ! IN299 TYPE(t_field), POINTER :: f_q(:) ! IN300 TYPE(t_field), POINTER :: f_temp(:) ! INOUT301 302 REAL(rstd), POINTER :: pmid(:,:)303 REAL(rstd), POINTER :: q(:,:,:)304 REAL(rstd), POINTER :: temp(:,:)305 INTEGER :: ind306 307 DO ind=1,ndomain308 IF (.NOT. assigned_domain(ind)) CYCLE309 CALL swap_dimensions(ind)310 CALL swap_geometry(ind)311 pmid=f_pmid(ind)312 q=f_q(ind)313 temp=f_temp(ind)314 CALL compute_diagnose_temp(pmid,q,temp)315 END DO316 END SUBROUTINE diagnose_temperature317 318 SUBROUTINE compute_diagnose_temp(pmid,q,temp)319 USE omp_para320 REAL(rstd),INTENT(IN) :: pmid(iim*jjm,llm)321 REAL(rstd),INTENT(IN) :: q(iim*jjm,llm,nqtot)322 REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm)323 324 REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix325 INTEGER :: ij,l326 327 Rd = kappa*cpp328 DO l=ll_begin,ll_end329 DO ij=ij_begin,ij_end330 p_ik = pmid(ij,l)331 theta_ik = temp(ij,l)332 qv = q(ij,l,1) ! water vapor mixing ratio = mv/md333 SELECT CASE(caldyn_thermo)334 CASE(thermo_theta)335 temp_ik = theta_ik*((p_ik/preff)**kappa)336 CASE(thermo_entropy)337 temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)338 CASE(thermo_moist)339 Rmix = Rd+qv*Rv340 chi = ( theta_ik + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)341 temp_ik = Treff*exp(chi)342 END SELECT343 IF(physics_thermo==thermo_fake_moist) temp_ik=temp_ik/(1+0.608*qv)344 temp(ij,l)=temp_ik345 END DO346 END DO347 END SUBROUTINE compute_diagnose_temp348 349 SUBROUTINE Tv2T(f_Tv, f_q, f_T)350 TYPE(t_field), POINTER :: f_TV(:)351 TYPE(t_field), POINTER :: f_q(:)352 TYPE(t_field), POINTER :: f_T(:)353 354 REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)355 INTEGER :: ind356 357 DO ind=1,ndomain358 IF (.NOT. assigned_domain(ind)) CYCLE359 CALL swap_dimensions(ind)360 CALL swap_geometry(ind)361 Tv=f_Tv(ind)362 T=f_T(ind)363 SELECT CASE(physics_thermo)364 CASE(thermo_dry)365 T=Tv366 CASE(thermo_fake_moist)367 q=f_q(ind)368 T=Tv/(1+0.608*q(:,:,1))369 END SELECT370 END DO371 END SUBROUTINE Tv2T372 373 214 SUBROUTINE divide_by_mass(iq, f_mass, f_theta_rhodz, f_theta) 374 215 INTEGER, INTENT(IN) :: iq
Note: See TracChangeset
for help on using the changeset viewer.