[325] | 1 | MODULE etat0_venus_mod |
---|
| 2 | USE icosa |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | PRIVATE |
---|
[741] | 5 | SAVE |
---|
[325] | 6 | |
---|
| 7 | TYPE(t_field),POINTER :: f_temp_eq( :) |
---|
| 8 | TYPE(t_field),POINTER :: f_temp(:) ! buffer used for physics |
---|
[741] | 9 | REAL(rstd), ALLOCATABLE :: temp_eq_packed(:,:) |
---|
[325] | 10 | |
---|
[741] | 11 | REAL(rstd) :: kfrict, kv |
---|
| 12 | !$OMP THREADPRIVATE(kfrict, kv) |
---|
[325] | 13 | |
---|
[741] | 14 | REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012 |
---|
| 15 | |
---|
| 16 | PUBLIC :: etat0, init_physics, full_physics |
---|
[325] | 17 | |
---|
| 18 | CONTAINS |
---|
| 19 | |
---|
| 20 | !-------------------------------- "Physics" ---------------------------------------- |
---|
| 21 | |
---|
[741] | 22 | SUBROUTINE init_physics_old |
---|
[325] | 23 | USE getin_mod |
---|
| 24 | REAL(rstd),POINTER :: temp(:,:) |
---|
| 25 | REAL(rstd) :: friction_time |
---|
| 26 | INTEGER :: ind |
---|
| 27 | |
---|
[741] | 28 | kv=0.15 ! vertical turbulent viscosity |
---|
| 29 | CALL getin('venus_diffusion',kv) |
---|
[325] | 30 | friction_time=86400. !friction |
---|
[741] | 31 | CALL getin('venus_friction_time',friction_time) |
---|
[325] | 32 | kfrict=1./friction_time |
---|
| 33 | |
---|
[741] | 34 | CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use by physics |
---|
[325] | 35 | |
---|
| 36 | PRINT *, 'Initializing Temp_eq (venus)' |
---|
| 37 | CALL allocate_field(f_temp_eq,field_t,type_real,llm) |
---|
| 38 | DO ind=1,ndomain |
---|
| 39 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
| 40 | CALL swap_dimensions(ind) |
---|
| 41 | CALL swap_geometry(ind) |
---|
| 42 | temp=f_temp_eq(ind) |
---|
[741] | 43 | CALL compute_temp_ref(.TRUE.,iim*jjm,lat_i, temp) ! With meridional gradient |
---|
[325] | 44 | ENDDO |
---|
[741] | 45 | END SUBROUTINE init_physics_old |
---|
[325] | 46 | |
---|
| 47 | SUBROUTINE physics(f_ps,f_theta_rhodz,f_u) |
---|
| 48 | USE theta2theta_rhodz_mod |
---|
| 49 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
---|
| 50 | TYPE(t_field),POINTER :: f_u(:) |
---|
| 51 | TYPE(t_field),POINTER :: f_ps(:) |
---|
| 52 | REAL(rstd),POINTER :: temp(:,:) |
---|
| 53 | REAL(rstd),POINTER :: temp_eq(:,:) |
---|
| 54 | REAL(rstd),POINTER :: u(:,:) |
---|
| 55 | INTEGER :: ind |
---|
| 56 | |
---|
| 57 | CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp) |
---|
| 58 | DO ind=1,ndomain |
---|
| 59 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
| 60 | CALL swap_dimensions(ind) |
---|
| 61 | CALL swap_geometry(ind) |
---|
| 62 | u=f_u(ind) |
---|
| 63 | temp_eq=f_temp_eq(ind) |
---|
| 64 | temp=f_temp(ind) |
---|
| 65 | CALL compute_physics(temp_eq, temp, u) |
---|
| 66 | ENDDO |
---|
| 67 | CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) |
---|
| 68 | END SUBROUTINE physics |
---|
| 69 | |
---|
| 70 | SUBROUTINE compute_physics(temp_eq, temp, u) |
---|
| 71 | REAL(rstd),INTENT(IN) :: temp_eq(iim*jjm,llm) |
---|
| 72 | REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm) |
---|
| 73 | REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm) |
---|
| 74 | INTEGER :: i,j,l,ij |
---|
| 75 | DO l=1,llm |
---|
| 76 | DO j=jj_begin-1,jj_end+1 |
---|
| 77 | DO i=ii_begin-1,ii_end+1 |
---|
| 78 | ij=(j-1)*iim+i |
---|
| 79 | temp(ij,l) = temp(ij,l) - (temp(ij,l)-temp_eq(ij,l))*(dt*itau_physics/tauCLee) |
---|
| 80 | ENDDO |
---|
| 81 | ENDDO |
---|
| 82 | ENDDO |
---|
| 83 | u(:,1)=u(:,1)*(1.-dt*itau_physics*kfrict) |
---|
| 84 | END SUBROUTINE compute_physics |
---|
| 85 | |
---|
[741] | 86 | !----------------------- Re-implementation using physics_interface_mod ----------------- |
---|
| 87 | |
---|
| 88 | SUBROUTINE init_physics |
---|
| 89 | USE getin_mod |
---|
| 90 | USE physics_interface_mod |
---|
| 91 | REAL(rstd),POINTER :: temp(:,:) |
---|
| 92 | REAL(rstd) :: friction_time |
---|
| 93 | INTEGER :: ngrid |
---|
| 94 | |
---|
| 95 | kv=0.15 ! vertical turbulent viscosity |
---|
| 96 | CALL getin('venus_diffusion',kv) |
---|
| 97 | friction_time=86400. !friction |
---|
| 98 | CALL getin('venus_friction_time',friction_time) |
---|
| 99 | kfrict=1./friction_time |
---|
| 100 | |
---|
| 101 | ngrid = physics_inout%ngrid |
---|
| 102 | ALLOCATE(temp_eq_packed(ngrid,llm)) |
---|
| 103 | CALL compute_temp_ref(.TRUE.,ngrid,physics_inout%lat, temp_eq_packed) ! With meridional gradient |
---|
| 104 | END SUBROUTINE init_physics |
---|
| 105 | |
---|
| 106 | SUBROUTINE full_physics |
---|
| 107 | USE physics_interface_mod |
---|
| 108 | CALL compute_physics_column(physics_inout%ngrid, physics_inout%dt_phys, & |
---|
| 109 | physics_inout%p, physics_inout%geopot, & |
---|
| 110 | physics_inout%Temp, physics_inout%ulon, physics_inout%ulat, & |
---|
| 111 | physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat) |
---|
| 112 | END SUBROUTINE full_physics |
---|
| 113 | |
---|
| 114 | SUBROUTINE compute_physics_column(ngrid,dt_phys,p,geopot,Temp,u,v,dTemp,du,dv) |
---|
| 115 | USE earth_const, only : g |
---|
| 116 | INTEGER, INTENT(IN) :: ngrid |
---|
| 117 | REAL(rstd),INTENT(IN) :: dt_phys, & |
---|
| 118 | p(ngrid,llm+1), geopot(ngrid,llm+1), & |
---|
| 119 | Temp(ngrid,llm), u(ngrid,llm), v(ngrid,llm) |
---|
| 120 | REAL(rstd),INTENT(OUT) :: dTemp(ngrid,llm), du(ngrid,llm), dv(ngrid,llm) |
---|
| 121 | REAL(rstd) :: mass(ngrid,llm), & ! rho.dz |
---|
| 122 | A(ngrid,llm+1), & ! off-diagonal coefficients |
---|
| 123 | B(ngrid,llm), & ! diagonal coefficients |
---|
| 124 | Ru(ngrid,llm), Rv(ngrid,llm), & ! right-hand-sides |
---|
| 125 | C(ngrid,llm), & ! LU factorization (Thomas algorithm) |
---|
| 126 | xu(ngrid,llm), xv(ngrid, llm) ! solution |
---|
| 127 | REAL(rstd) :: rho, X_ij |
---|
| 128 | INTEGER :: l,ij |
---|
| 129 | ! Vertical diffusion : rho.du/dt = d/dz (rho*kappa*du/dz) |
---|
| 130 | ! rho.dz = mass.deta |
---|
| 131 | ! => mass.du/dt = d/deta ((rho^2 kappa/m)du/deta) |
---|
| 132 | ! Backward Euler : (M-S)u_new = M.u_old |
---|
| 133 | ! with M.u = mass(l)*u(l) |
---|
| 134 | ! and S.u = A(l+1)*(u(l+1)-u(l)) - A(l)*(u(l)-u(l-1)) |
---|
| 135 | ! => solve -A(l)u(l-1) + B(l)u(l) - A(l+1)u(l+1) = Ru(l) using Thomas algorithm |
---|
| 136 | ! with A=tau*kappa*rho^2/m and B(l) = mass(l)+A(l)+A(l+1) |
---|
| 137 | DO l=1,llm |
---|
| 138 | !DIR$ SIMD |
---|
| 139 | DO ij=1,ngrid |
---|
| 140 | mass(ij,l) = (p(ij,l)-p(ij,l+1))*(1./g) |
---|
| 141 | rho = (p(ij,l)-p(ij,l+1))/(geopot(ij,l+1)-geopot(ij,l)) |
---|
| 142 | ! A = kappa.tau.rho^2/m |
---|
| 143 | A(ij,l) = (kv*dt_phys)*(rho**2)/mass(ij,l) |
---|
| 144 | Ru(ij,l) = mass(ij,l)*u(ij,l) |
---|
| 145 | Rv(ij,l) = mass(ij,l)*v(ij,l) |
---|
| 146 | ENDDO |
---|
| 147 | ENDDO |
---|
| 148 | A(:,llm+1)=0. |
---|
| 149 | DO l=llm,2,-1 |
---|
| 150 | !DIR$ SIMD |
---|
| 151 | DO ij=1,ngrid |
---|
| 152 | A(ij,l) = .5*(A(ij,l)+A(ij,l-1)) ! average A to interfaces |
---|
| 153 | ENDDO |
---|
| 154 | ENDDO |
---|
| 155 | A(:,1)=0. |
---|
| 156 | DO l=1,llm |
---|
| 157 | !DIR$ SIMD |
---|
| 158 | DO ij=1,ngrid |
---|
| 159 | B(ij,l) = mass(ij,l)+A(ij,l)+A(ij,l+1) |
---|
| 160 | ENDDO |
---|
| 161 | ENDDO |
---|
| 162 | |
---|
| 163 | ! Solve -A(l)x(l-1) + B(l)x(l) - A(l+1)x(l+1) = R(l) using Thomas algorithm |
---|
| 164 | ! Forward sweep : |
---|
| 165 | ! C(0)=0, C(l) = -A(l+1) / (B(l)+A(l)C(l-1)), |
---|
| 166 | ! D(0)=0, D(l) = (R(l)+A(l)D(l-1)) / (B(l)+A(l)C(l-1)) |
---|
| 167 | !DIR$ SIMD |
---|
| 168 | DO ij=1,ngrid |
---|
| 169 | X_ij = 1./B(ij,1) |
---|
| 170 | C(ij,2) = -A(ij,2) * X_ij |
---|
| 171 | xu(ij,1) = Ru(ij,1) * X_ij |
---|
| 172 | xv(ij,1) = Rv(ij,1) * X_ij |
---|
| 173 | END DO |
---|
| 174 | DO l = 2,llm |
---|
| 175 | !DIR$ SIMD |
---|
| 176 | DO ij=1,ngrid |
---|
| 177 | X_ij = 1./( B(ij,l) + A(ij,l)*C(ij,l) ) |
---|
| 178 | C(ij,l+1) = -A(ij,l+1) * X_ij ! zero for l=llm |
---|
| 179 | xu(ij,l) = (Ru(ij,l)+A(ij,l)*xu(ij,l-1)) * X_ij |
---|
| 180 | xv(ij,l) = (Rv(ij,l)+A(ij,l)*xv(ij,l-1)) * X_ij |
---|
| 181 | END DO |
---|
| 182 | END DO |
---|
| 183 | ! Back substitution : |
---|
| 184 | ! x(i) = D(i)-C(i+1)x(i+1), x(llm)=D(llm) |
---|
| 185 | !DIR$ SIMD |
---|
| 186 | DO ij=1,ngrid |
---|
| 187 | ! top layer l=llm |
---|
| 188 | du(ij,llm) = (xu(ij,llm)-u(ij,llm))/dt_phys |
---|
| 189 | dv(ij,llm) = (xv(ij,llm)-v(ij,llm))/dt_phys |
---|
| 190 | END DO |
---|
| 191 | ! Back substitution at lower layers |
---|
| 192 | DO l = llm-1,1,-1 |
---|
| 193 | !DIR$ SIMD |
---|
| 194 | DO ij=1,ngrid |
---|
| 195 | xu(ij,l) = xu(ij,l) - C(ij,l+1)*xu(ij,l+1) |
---|
| 196 | xv(ij,l) = xv(ij,l) - C(ij,l+1)*xv(ij,l+1) |
---|
| 197 | du(ij,l) = (xu(ij,l)-u(ij,l))/dt_phys |
---|
| 198 | dv(ij,l) = (xv(ij,l)-v(ij,l))/dt_phys |
---|
| 199 | END DO |
---|
| 200 | END DO |
---|
| 201 | ! bottom friction + thermal relaxation |
---|
| 202 | du(:,1)=du(:,1)-kfrict*u(:,1) |
---|
| 203 | dTemp(:,:) = (1./tauCLee)*(temp_eq_packed(:,:)-temp(:,:)) |
---|
| 204 | END SUBROUTINE compute_physics_column |
---|
| 205 | |
---|
[325] | 206 | !----------------------------- Initialize to T_eq -------------------------------------- |
---|
| 207 | |
---|
| 208 | SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) |
---|
| 209 | USE theta2theta_rhodz_mod |
---|
| 210 | TYPE(t_field),POINTER :: f_ps(:) |
---|
| 211 | TYPE(t_field),POINTER :: f_phis(:) |
---|
| 212 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
---|
| 213 | TYPE(t_field),POINTER :: f_u(:) |
---|
| 214 | TYPE(t_field),POINTER :: f_q(:) |
---|
| 215 | |
---|
| 216 | TYPE(t_field),POINTER :: f_temp(:) |
---|
| 217 | REAL(rstd),POINTER :: temp(:,:) |
---|
| 218 | REAL(rstd),POINTER :: ps(:) |
---|
| 219 | REAL(rstd),POINTER :: phis(:) |
---|
| 220 | REAL(rstd),POINTER :: u(:,:) |
---|
| 221 | REAL(rstd),POINTER :: q(:,:,:) |
---|
| 222 | INTEGER :: ind |
---|
| 223 | |
---|
[741] | 224 | CALL allocate_field(f_temp,field_t,type_real,llm, name='temp_buf_venus') |
---|
[325] | 225 | DO ind=1,ndomain |
---|
| 226 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
| 227 | CALL swap_dimensions(ind) |
---|
| 228 | CALL swap_geometry(ind) |
---|
| 229 | ps=f_ps(ind) |
---|
| 230 | ps(:)=preff |
---|
| 231 | phis=f_phis(ind) |
---|
| 232 | phis(:)=0. |
---|
| 233 | u=f_u(ind) |
---|
| 234 | u(:,:)=0 |
---|
| 235 | q=f_q(ind) |
---|
| 236 | q(:,:,:)=1e2 |
---|
| 237 | temp=f_temp(ind) |
---|
[741] | 238 | CALL compute_temp_ref(.FALSE., iim*jjm, lat_i, temp) ! Without meridional gradient |
---|
[325] | 239 | ENDDO |
---|
| 240 | CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) |
---|
| 241 | CALL deallocate_field(f_temp) |
---|
| 242 | |
---|
| 243 | END SUBROUTINE etat0 |
---|
| 244 | |
---|
| 245 | !------------------------- Compute reference temperature field ------------------------ |
---|
| 246 | |
---|
[741] | 247 | SUBROUTINE compute_temp_ref(gradient,ngrid,lat, temp_eq) |
---|
[325] | 248 | USE disvert_mod |
---|
| 249 | USE omp_para |
---|
| 250 | USE math_const |
---|
| 251 | LOGICAL, INTENT(IN) :: gradient |
---|
[741] | 252 | INTEGER, INTENT(IN) :: ngrid |
---|
| 253 | REAL(rstd), INTENT(IN) :: lat(ngrid) ! latitude |
---|
| 254 | REAL(rstd), INTENT(OUT) :: temp_eq(ngrid,llm) |
---|
| 255 | |
---|
| 256 | REAL(rstd) :: clat(ngrid) |
---|
| 257 | INTEGER :: level |
---|
| 258 | INTEGER, PARAMETER :: nlev=30 |
---|
| 259 | REAL :: pressCLee(nlev+1), tempCLee(nlev+1), dt_epCLee(nlev+1), etaCLee(nlev+1) |
---|
[325] | 260 | |
---|
[741] | 261 | REAL(rstd) :: pplay, ztemp,zdt,fact |
---|
| 262 | INTEGER :: ij, l,ll |
---|
[325] | 263 | |
---|
| 264 | data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, & |
---|
| 265 | 4.273e-1, 3.373e-1, 2.610e-1,1.979e-1,1.472e-1, & |
---|
| 266 | 1.074e-1, 7.672e-2, 5.361e-2,3.657e-2,2.430e-2, & |
---|
| 267 | 1.569e-2, 9.814e-3, 5.929e-3,3.454e-3,1.934e-3, & |
---|
| 268 | 1.043e-3, 5.400e-4, 2.710e-4,1.324e-4,6.355e-5, & |
---|
| 269 | 3.070e-5, 1.525e-5, 7.950e-6,4.500e-6,2.925e-6, & |
---|
| 270 | 2.265e-6/ |
---|
| 271 | |
---|
| 272 | |
---|
| 273 | data tempCLee/ 728.187, 715.129, 697.876, 677.284, 654.078, 628.885, & |
---|
| 274 | 602.225, 574.542, 546.104, 517.339, 488.560, 459.932, & |
---|
| 275 | 431.741, 404.202, 377.555, 352.042, 327.887, 305.313, & |
---|
| 276 | 284.556, 265.697, 248.844, 233.771, 220.368, 208.247, & |
---|
| 277 | 197.127, 187.104, 178.489, 171.800, 167.598, 165.899, & |
---|
| 278 | 165.676/ |
---|
| 279 | data dt_epCLee/6.101 , 6.136 , 6.176 , 6.410 , 6.634 , 6.678 , & |
---|
| 280 | 6.719 , 6.762 , 7.167 , 7.524 , 9.840 ,14.948 , & |
---|
| 281 | 21.370 ,28.746 ,36.373 ,43.315 ,48.534 ,51.175 , & |
---|
| 282 | 50.757 ,47.342 ,41.536 ,34.295 ,26.758 ,19.807 , & |
---|
| 283 | 14.001 , 9.599 , 6.504 , 4.439 , 3.126 , 2.370 , & |
---|
| 284 | 2.000/ |
---|
| 285 | |
---|
[741] | 286 | pressCLee(:) = etaCLee(:)*9.2e6 |
---|
| 287 | clat(:)=COS(lat(:)) |
---|
[325] | 288 | |
---|
[741] | 289 | DO ij=1,ngrid |
---|
| 290 | DO l = 1, llm |
---|
| 291 | |
---|
| 292 | pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff |
---|
| 293 | ! look for largest level such that pressCLee(level) > pplay(ij,l)) |
---|
| 294 | ! => pressClee(level+1) < pplay(ij,l) < pressClee(level) |
---|
| 295 | level = 1 |
---|
| 296 | DO ll = 1, nlev ! nlev data levels |
---|
| 297 | IF(pressCLee(ll) > pplay) THEN |
---|
| 298 | level = ll |
---|
[325] | 299 | END IF |
---|
| 300 | END DO |
---|
[741] | 301 | |
---|
| 302 | ! interpolate between level and level+1 |
---|
| 303 | ! interpolation is linear in log(pressure) |
---|
| 304 | fact = ( log10(pplay)-log10(pressCLee(level))) & |
---|
| 305 | /( log10(pressCLee(level+1))-log10(pressCLee(level)) ) |
---|
| 306 | ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact |
---|
| 307 | zdt = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact |
---|
| 308 | |
---|
| 309 | IF(gradient) THEN |
---|
| 310 | temp_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.) |
---|
| 311 | ELSE |
---|
| 312 | temp_eq(ij,l) = ztemp |
---|
| 313 | END IF |
---|
[325] | 314 | END DO |
---|
| 315 | END DO |
---|
| 316 | END SUBROUTINE compute_temp_ref |
---|
| 317 | |
---|
| 318 | END MODULE etat0_venus_mod |
---|