1 | MODULE etat0_venus_mod |
---|
2 | USE icosa |
---|
3 | IMPLICIT NONE |
---|
4 | PRIVATE |
---|
5 | SAVE |
---|
6 | |
---|
7 | TYPE(t_field),POINTER :: f_temp_eq( :) |
---|
8 | TYPE(t_field),POINTER :: f_temp(:) ! buffer used for physics |
---|
9 | REAL(rstd), ALLOCATABLE :: temp_eq_packed(:,:) |
---|
10 | |
---|
11 | REAL(rstd) :: kfrict, kv |
---|
12 | !$OMP THREADPRIVATE(kfrict, kv) |
---|
13 | |
---|
14 | REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012 |
---|
15 | |
---|
16 | PUBLIC :: etat0, init_physics, full_physics |
---|
17 | |
---|
18 | CONTAINS |
---|
19 | |
---|
20 | !-------------------------------- "Physics" ---------------------------------------- |
---|
21 | |
---|
22 | SUBROUTINE init_physics_old |
---|
23 | USE getin_mod |
---|
24 | REAL(rstd),POINTER :: temp(:,:) |
---|
25 | REAL(rstd) :: friction_time |
---|
26 | INTEGER :: ind |
---|
27 | |
---|
28 | kv=0.15 ! vertical turbulent viscosity |
---|
29 | CALL getin('venus_diffusion',kv) |
---|
30 | friction_time=86400. !friction |
---|
31 | CALL getin('venus_friction_time',friction_time) |
---|
32 | kfrict=1./friction_time |
---|
33 | |
---|
34 | CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use by physics |
---|
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) |
---|
43 | CALL compute_temp_ref(.TRUE.,iim*jjm,lat_i, temp) ! With meridional gradient |
---|
44 | ENDDO |
---|
45 | END SUBROUTINE init_physics_old |
---|
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 | |
---|
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,1) = -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-1) ) |
---|
178 | C(ij,l) = -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)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)*xu(ij,l+1) |
---|
196 | xv(ij,l) = xv(ij,l) - C(ij,l)*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 | |
---|
202 | IF(.FALSE.) THEN |
---|
203 | ! check correctness of solution of tridiagonal problem |
---|
204 | PRINT *, 'Max 1 residual of tridiagonal solver', MAXVAL(ABS(Ru)) |
---|
205 | DO l = 2,llm-1 |
---|
206 | !DIR$ SIMD |
---|
207 | DO ij=1,ngrid |
---|
208 | Ru(ij,l) = Ru(ij,l) - B(ij,l)*xu(ij,l) + A(ij,l)*xu(ij,l-1) + A(ij,l+1)*xu(ij,l+1) |
---|
209 | END DO |
---|
210 | END DO |
---|
211 | DO ij=1,ngrid |
---|
212 | Ru(ij,1) = Ru(ij,1) - B(ij,1)*xu(ij,1) + A(ij,2)*xu(ij,2) |
---|
213 | Ru(ij,llm) = Ru(ij,llm) - B(ij,llm)*xu(ij,llm) + A(ij,llm)*xu(ij,llm-1) |
---|
214 | END DO |
---|
215 | PRINT *, 'Max 2 residual of tridiagonal solver', MAXVAL(ABS(Ru)) |
---|
216 | END IF |
---|
217 | |
---|
218 | ! bottom friction + thermal relaxation |
---|
219 | du(:,1)=du(:,1)-kfrict*u(:,1) |
---|
220 | dTemp(:,:) = (1./tauCLee)*(temp_eq_packed(:,:)-temp(:,:)) |
---|
221 | END SUBROUTINE compute_physics_column |
---|
222 | |
---|
223 | !----------------------------- Initialize to T_eq -------------------------------------- |
---|
224 | |
---|
225 | SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) |
---|
226 | USE theta2theta_rhodz_mod |
---|
227 | TYPE(t_field),POINTER :: f_ps(:) |
---|
228 | TYPE(t_field),POINTER :: f_phis(:) |
---|
229 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
---|
230 | TYPE(t_field),POINTER :: f_u(:) |
---|
231 | TYPE(t_field),POINTER :: f_q(:) |
---|
232 | |
---|
233 | TYPE(t_field),POINTER :: f_temp(:) |
---|
234 | REAL(rstd),POINTER :: temp(:,:) |
---|
235 | REAL(rstd),POINTER :: ps(:) |
---|
236 | REAL(rstd),POINTER :: phis(:) |
---|
237 | REAL(rstd),POINTER :: u(:,:) |
---|
238 | REAL(rstd),POINTER :: q(:,:,:) |
---|
239 | INTEGER :: ind |
---|
240 | |
---|
241 | CALL allocate_field(f_temp,field_t,type_real,llm, name='temp_buf_venus') |
---|
242 | DO ind=1,ndomain |
---|
243 | IF (.NOT. assigned_domain(ind)) CYCLE |
---|
244 | CALL swap_dimensions(ind) |
---|
245 | CALL swap_geometry(ind) |
---|
246 | ps=f_ps(ind) |
---|
247 | ps(:)=preff |
---|
248 | phis=f_phis(ind) |
---|
249 | phis(:)=0. |
---|
250 | u=f_u(ind) |
---|
251 | u(:,:)=0 |
---|
252 | q=f_q(ind) |
---|
253 | q(:,:,:)=1e2 |
---|
254 | temp=f_temp(ind) |
---|
255 | CALL compute_temp_ref(.FALSE., iim*jjm, lat_i, temp) ! Without meridional gradient |
---|
256 | ENDDO |
---|
257 | CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) |
---|
258 | CALL deallocate_field(f_temp) |
---|
259 | |
---|
260 | END SUBROUTINE etat0 |
---|
261 | |
---|
262 | !------------------------- Compute reference temperature field ------------------------ |
---|
263 | |
---|
264 | SUBROUTINE compute_temp_ref(gradient,ngrid,lat, temp_eq) |
---|
265 | USE disvert_mod |
---|
266 | USE omp_para |
---|
267 | USE math_const |
---|
268 | LOGICAL, INTENT(IN) :: gradient |
---|
269 | INTEGER, INTENT(IN) :: ngrid |
---|
270 | REAL(rstd), INTENT(IN) :: lat(ngrid) ! latitude |
---|
271 | REAL(rstd), INTENT(OUT) :: temp_eq(ngrid,llm) |
---|
272 | |
---|
273 | REAL(rstd) :: clat(ngrid) |
---|
274 | INTEGER :: level |
---|
275 | INTEGER, PARAMETER :: nlev=30 |
---|
276 | REAL :: pressCLee(nlev+1), tempCLee(nlev+1), dt_epCLee(nlev+1), etaCLee(nlev+1) |
---|
277 | |
---|
278 | REAL(rstd) :: pplay, ztemp,zdt,fact |
---|
279 | INTEGER :: ij, l,ll |
---|
280 | |
---|
281 | data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, & |
---|
282 | 4.273e-1, 3.373e-1, 2.610e-1,1.979e-1,1.472e-1, & |
---|
283 | 1.074e-1, 7.672e-2, 5.361e-2,3.657e-2,2.430e-2, & |
---|
284 | 1.569e-2, 9.814e-3, 5.929e-3,3.454e-3,1.934e-3, & |
---|
285 | 1.043e-3, 5.400e-4, 2.710e-4,1.324e-4,6.355e-5, & |
---|
286 | 3.070e-5, 1.525e-5, 7.950e-6,4.500e-6,2.925e-6, & |
---|
287 | 2.265e-6/ |
---|
288 | |
---|
289 | |
---|
290 | data tempCLee/ 728.187, 715.129, 697.876, 677.284, 654.078, 628.885, & |
---|
291 | 602.225, 574.542, 546.104, 517.339, 488.560, 459.932, & |
---|
292 | 431.741, 404.202, 377.555, 352.042, 327.887, 305.313, & |
---|
293 | 284.556, 265.697, 248.844, 233.771, 220.368, 208.247, & |
---|
294 | 197.127, 187.104, 178.489, 171.800, 167.598, 165.899, & |
---|
295 | 165.676/ |
---|
296 | data dt_epCLee/6.101 , 6.136 , 6.176 , 6.410 , 6.634 , 6.678 , & |
---|
297 | 6.719 , 6.762 , 7.167 , 7.524 , 9.840 ,14.948 , & |
---|
298 | 21.370 ,28.746 ,36.373 ,43.315 ,48.534 ,51.175 , & |
---|
299 | 50.757 ,47.342 ,41.536 ,34.295 ,26.758 ,19.807 , & |
---|
300 | 14.001 , 9.599 , 6.504 , 4.439 , 3.126 , 2.370 , & |
---|
301 | 2.000/ |
---|
302 | |
---|
303 | pressCLee(:) = etaCLee(:)*9.2e6 |
---|
304 | clat(:)=COS(lat(:)) |
---|
305 | |
---|
306 | DO ij=1,ngrid |
---|
307 | DO l = 1, llm |
---|
308 | |
---|
309 | pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff |
---|
310 | ! look for largest level such that pressCLee(level) > pplay(ij,l)) |
---|
311 | ! => pressClee(level+1) < pplay(ij,l) < pressClee(level) |
---|
312 | level = 1 |
---|
313 | DO ll = 1, nlev ! nlev data levels |
---|
314 | IF(pressCLee(ll) > pplay) THEN |
---|
315 | level = ll |
---|
316 | END IF |
---|
317 | END DO |
---|
318 | |
---|
319 | ! interpolate between level and level+1 |
---|
320 | ! interpolation is linear in log(pressure) |
---|
321 | fact = ( log10(pplay)-log10(pressCLee(level))) & |
---|
322 | /( log10(pressCLee(level+1))-log10(pressCLee(level)) ) |
---|
323 | ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact |
---|
324 | zdt = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact |
---|
325 | |
---|
326 | IF(gradient) THEN |
---|
327 | temp_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.) |
---|
328 | ELSE |
---|
329 | temp_eq(ij,l) = ztemp |
---|
330 | END IF |
---|
331 | END DO |
---|
332 | END DO |
---|
333 | END SUBROUTINE compute_temp_ref |
---|
334 | |
---|
335 | END MODULE etat0_venus_mod |
---|