[377] | 1 | MODULE etat0_bubble_mod |
---|
| 2 | USE icosa |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | PRIVATE |
---|
| 5 | |
---|
| 6 | REAL(rstd), PARAMETER :: latc=0., lonc=0., & |
---|
| 7 | width = 0.3, & ! width of mountain (radian) |
---|
| 8 | width_bubble = 0.5, & ! width of bubble (radian), curvature effects O(width^2) |
---|
| 9 | zc = 350., & ! height at bubble center |
---|
| 10 | zw = 250. ! bubble vertical scale => radius=4000m for circular bubble |
---|
| 11 | REAL(rstd), SAVE :: T0 ! adiabatic atmosphere at theta=T0 (default 300K) |
---|
| 12 | !$OMP THREADPRIVATE(T0) |
---|
| 13 | REAL(rstd), SAVE :: ztop ! mountain top (default 0.) |
---|
| 14 | !$OMP THREADPRIVATE(ztop) |
---|
| 15 | REAL(rstd), SAVE :: delta_T ! bubble extra temperature (default 5K) |
---|
| 16 | !$OMP THREADPRIVATE(delta_T) |
---|
| 17 | |
---|
| 18 | PUBLIC getin_etat0, compute_etat0 |
---|
| 19 | |
---|
| 20 | CONTAINS |
---|
| 21 | |
---|
| 22 | SUBROUTINE getin_etat0 |
---|
| 23 | T0=300. |
---|
| 24 | CALL getin("bubble_T0",T0) |
---|
| 25 | ztop=0. |
---|
| 26 | CALL getin("bubble_mountain_top",ztop) |
---|
| 27 | delta_T=5. |
---|
| 28 | CALL getin("bubble_delta_T",delta_T) |
---|
| 29 | END SUBROUTINE getin_etat0 |
---|
| 30 | |
---|
| 31 | SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,geopot,q) |
---|
| 32 | USE icosa |
---|
| 33 | USE disvert_mod |
---|
| 34 | USE omp_para |
---|
| 35 | INTEGER, INTENT(IN) :: ngrid |
---|
| 36 | REAL(rstd),INTENT(IN) :: lon(ngrid) |
---|
| 37 | REAL(rstd),INTENT(IN) :: lat(ngrid) |
---|
| 38 | REAL(rstd),INTENT(OUT) :: phis(ngrid) |
---|
| 39 | REAL(rstd),INTENT(OUT) :: ps(ngrid) |
---|
| 40 | REAL(rstd),INTENT(OUT) :: temp(ngrid,llm) |
---|
| 41 | REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) |
---|
| 42 | REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) |
---|
| 43 | REAL(rstd),INTENT(OUT) :: geopot(ngrid,llm+1) |
---|
| 44 | REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) |
---|
| 45 | |
---|
| 46 | INTEGER :: l,ij |
---|
| 47 | REAL(rstd) :: sinlat, coslat, K, Ts, p_ij, T_ij, z_ij |
---|
| 48 | ! adiabatic profile |
---|
| 49 | ! p=preff, T=T0 at Phi=0 |
---|
| 50 | ! => theta = T0, cp.T + Phi = cst = cp.T0 |
---|
| 51 | DO ij=1,ngrid |
---|
| 52 | sinlat=SIN(lat(ij)) |
---|
| 53 | coslat=COS(lat(ij)) |
---|
| 54 | K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc) |
---|
| 55 | K=acos(K)/width |
---|
| 56 | phis(ij) = g*ztop*exp(-(K**2)) |
---|
| 57 | Ts = T0-phis(ij)/cpp ! adiabatic |
---|
| 58 | ps(ij) = preff * (Ts/T0)**(1/kappa) ! theta=T0=T.(p/preff)^-kappa |
---|
| 59 | ENDDO |
---|
| 60 | |
---|
| 61 | DO l=1,llm+1 |
---|
| 62 | DO ij=1,ngrid |
---|
| 63 | p_ij = ap(l)+bp(l)*ps(ij) |
---|
| 64 | T_ij = T0*(p_ij/preff)**kappa ! adiabatic |
---|
| 65 | geopot(ij,l) = cpp*(T0-T_ij) |
---|
| 66 | END DO |
---|
| 67 | END DO |
---|
| 68 | |
---|
| 69 | DO l=1, llm |
---|
| 70 | DO ij=1,ngrid |
---|
| 71 | sinlat=SIN(lat(ij)) |
---|
| 72 | coslat=COS(lat(ij)) |
---|
| 73 | K=sin(latc)*sinlat+cos(latc)*coslat |
---|
| 74 | K=acos(K)/width_bubble |
---|
| 75 | p_ij = 0.5 *( ap(l)+bp(l)*ps(ij) + ap(l+1)+bp(l+1)*ps(ij) ) |
---|
| 76 | T_ij = T0*(p_ij/preff)**kappa |
---|
| 77 | z_ij = (T0-T_ij)*cpp/g |
---|
| 78 | z_ij = (z_ij - zc)/zw |
---|
| 79 | K=MIN(1.,sqrt(K**2+z_ij**2)) |
---|
| 80 | K=.5*(1.+COS(pi*K)) |
---|
| 81 | Temp(ij,l) = T_ij + delta_T*K |
---|
| 82 | ulat(ij,l)=0. |
---|
| 83 | ulon(ij,l)=0. |
---|
| 84 | q(ij,l,:)=0. |
---|
| 85 | END DO |
---|
| 86 | END DO |
---|
| 87 | |
---|
| 88 | END SUBROUTINE compute_etat0 |
---|
| 89 | |
---|
| 90 | END MODULE etat0_bubble_mod |
---|