source: codes/icosagcm/devel/src/initial/etat0_venus.f90 @ 847

Last change on this file since 847 was 741, checked in by dubos, 6 years ago

devel : added vertical diffusion to idealized venus physics
File size: 10.9 KB
Line 
1MODULE 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 
18CONTAINS
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,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
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
224    CALL allocate_field(f_temp,field_t,type_real,llm, name='temp_buf_venus')
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)
238       CALL compute_temp_ref(.FALSE., iim*jjm, lat_i, temp) ! Without meridional gradient
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
247  SUBROUTINE compute_temp_ref(gradient,ngrid,lat, temp_eq)
248    USE disvert_mod
249    USE omp_para
250    USE math_const
251    LOGICAL, INTENT(IN) :: gradient
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)
260   
261    REAL(rstd) ::  pplay, ztemp,zdt,fact
262    INTEGER :: ij, l,ll
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     
286    pressCLee(:) = etaCLee(:)*9.2e6
287    clat(:)=COS(lat(:))
288
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
299             END IF
300          END DO
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
314       END DO
315    END DO
316  END SUBROUTINE compute_temp_ref
317
318END MODULE etat0_venus_mod
Note: See TracBrowser for help on using the repository browser.