source: codes/icosagcm/trunk/src/initial/etat0_venus.f90 @ 899

Last change on this file since 899 was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

Fixed iso c bindings
fixed warnings with -Wall -Wno-aliasing -Wno-unused -Wno-unused-dummy-argument -Wno-maybe-uninitialized -Wno-tabs warnings
Removed all unused variables (-Wunused-variable)
vector%dot_product is now dot_product_3d to avoid compilation warning "dot_product shadows intrinsic" with GCC
  • Property svn:executable set to *
File size: 6.8 KB
Line 
1MODULE etat0_venus_mod
2  USE icosa
3  IMPLICIT NONE
4  PRIVATE
5
6  TYPE(t_field),POINTER :: f_temp_eq( :)
7  TYPE(t_field),POINTER :: f_temp(:) ! buffer used for physics
8
9  REAL(rstd), SAVE :: kfrict
10!$OMP THREADPRIVATE(kfrict)
11
12  PUBLIC :: etat0, init_physics, physics
13 
14CONTAINS
15
16!-------------------------------- "Physics" ----------------------------------------
17
18  SUBROUTINE init_physics
19    USE getin_mod
20    IMPLICIT NONE
21    REAL(rstd),POINTER :: temp(:,:)
22    REAL(rstd) :: friction_time
23    INTEGER :: ind
24
25    friction_time=86400.  !friction
26    CALL getin('friction_time',friction_time)
27    kfrict=1./friction_time
28
29    CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use
30
31    PRINT *, 'Initializing Temp_eq (venus)'
32    CALL allocate_field(f_temp_eq,field_t,type_real,llm)
33    DO ind=1,ndomain
34       IF (.NOT. assigned_domain(ind)) CYCLE
35       CALL swap_dimensions(ind)
36       CALL swap_geometry(ind)
37       temp=f_temp_eq(ind)
38       CALL compute_temp_ref(temp, .TRUE.) ! FIXME With meridional gradient
39    ENDDO
40  END SUBROUTINE init_physics
41
42  SUBROUTINE physics(f_ps,f_theta_rhodz,f_u) 
43    USE icosa
44    USE theta2theta_rhodz_mod
45    IMPLICIT NONE
46    TYPE(t_field),POINTER :: f_theta_rhodz(:)
47    TYPE(t_field),POINTER :: f_u(:)
48    TYPE(t_field),POINTER :: f_ps(:)
49    REAL(rstd),POINTER :: temp(:,:)
50    REAL(rstd),POINTER :: temp_eq(:,:)
51    REAL(rstd),POINTER :: u(:,:)
52    INTEGER :: ind
53
54    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
55    DO ind=1,ndomain
56       IF (.NOT. assigned_domain(ind)) CYCLE
57       CALL swap_dimensions(ind)
58       CALL swap_geometry(ind)
59       u=f_u(ind)
60       temp_eq=f_temp_eq(ind) 
61       temp=f_temp(ind) 
62       CALL compute_physics(temp_eq, temp, u)
63    ENDDO
64    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
65  END SUBROUTINE physics
66
67  SUBROUTINE compute_physics(temp_eq, temp, u)
68    USE icosa
69    USE theta2theta_rhodz_mod
70    IMPLICIT NONE
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    REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012
75    INTEGER :: i,j,l,ij
76    DO l=1,llm
77       DO j=jj_begin-1,jj_end+1
78          DO i=ii_begin-1,ii_end+1
79             ij=(j-1)*iim+i
80             temp(ij,l) = temp(ij,l) - (temp(ij,l)-temp_eq(ij,l))*(dt*itau_physics/tauCLee)
81          ENDDO
82       ENDDO
83    ENDDO
84    u(:,1)=u(:,1)*(1.-dt*itau_physics*kfrict)
85  END SUBROUTINE compute_physics
86
87!----------------------------- Initialize to T_eq --------------------------------------
88
89  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
90    USE icosa
91    USE theta2theta_rhodz_mod
92    IMPLICIT NONE
93    TYPE(t_field),POINTER :: f_ps(:)
94    TYPE(t_field),POINTER :: f_phis(:)
95    TYPE(t_field),POINTER :: f_theta_rhodz(:)
96    TYPE(t_field),POINTER :: f_u(:)
97    TYPE(t_field),POINTER :: f_q(:)
98
99    TYPE(t_field),POINTER :: f_temp(:)
100    REAL(rstd),POINTER :: temp(:,:) 
101    REAL(rstd),POINTER :: ps(:)
102    REAL(rstd),POINTER :: phis(:)
103    REAL(rstd),POINTER :: u(:,:)
104    REAL(rstd),POINTER :: q(:,:,:)                 
105    INTEGER :: ind
106
107    CALL allocate_field(f_temp,field_t,type_real,llm)
108    DO ind=1,ndomain
109       IF (.NOT. assigned_domain(ind)) CYCLE
110       CALL swap_dimensions(ind)
111       CALL swap_geometry(ind)
112       ps=f_ps(ind)
113       ps(:)=preff 
114       phis=f_phis(ind)
115       phis(:)=0.
116       u=f_u(ind)
117       u(:,:)=0
118       q=f_q(ind)
119       q(:,:,:)=1e2
120       temp=f_temp(ind)
121       CALL compute_temp_ref(temp, .FALSE.) ! Without meridional gradient
122    ENDDO
123    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
124    CALL deallocate_field(f_temp)
125
126  END SUBROUTINE etat0
127
128!------------------------- Compute reference temperature field ------------------------
129
130  SUBROUTINE compute_temp_ref(theta_eq, gradient)
131    USE icosa
132    USE disvert_mod
133    USE omp_para
134    USE math_const
135    IMPLICIT NONE
136    REAL(rstd), INTENT(OUT) :: theta_eq(iim*jjm,llm) 
137    LOGICAL, INTENT(IN) :: gradient
138    REAL(rstd)  :: clat(iim*jjm)        ! latitude
139    integer :: level
140    REAL ::  pressCLee(31), tempCLee(31), dt_epCLee(31), etaCLee(31)
141   
142    real(rstd) ::  lon,lat, pplay, ztemp,zdt,fact
143    integer :: i,j,ij, l,ll
144   
145    data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, &
146                      4.273e-1, 3.373e-1, 2.610e-1,1.979e-1,1.472e-1,  &
147                      1.074e-1, 7.672e-2, 5.361e-2,3.657e-2,2.430e-2,  &
148                      1.569e-2, 9.814e-3, 5.929e-3,3.454e-3,1.934e-3,  &
149                      1.043e-3, 5.400e-4, 2.710e-4,1.324e-4,6.355e-5,  &
150                      3.070e-5, 1.525e-5, 7.950e-6,4.500e-6,2.925e-6,  &
151                      2.265e-6/
152
153    DO j=jj_begin-1,jj_end+1
154       DO i=ii_begin-1,ii_end+1
155          ij=(j-1)*iim+i
156          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
157          clat(ij)=cos(lat) 
158       ENDDO
159    ENDDO
160
161    data   tempCLee/ 728.187, 715.129, 697.876, 677.284, 654.078, 628.885, &
162                     602.225, 574.542, 546.104, 517.339, 488.560, 459.932, &
163                     431.741, 404.202, 377.555, 352.042, 327.887, 305.313, &
164                     284.556, 265.697, 248.844, 233.771, 220.368, 208.247, &
165                     197.127, 187.104, 178.489, 171.800, 167.598, 165.899, &
166                     165.676/
167    data   dt_epCLee/6.101 , 6.136 , 6.176 , 6.410 , 6.634 , 6.678 , &
168                     6.719 , 6.762 , 7.167 , 7.524 , 9.840 ,14.948 ,  &
169                     21.370 ,28.746 ,36.373 ,43.315 ,48.534 ,51.175 ,  &
170                     50.757 ,47.342 ,41.536 ,34.295 ,26.758 ,19.807 ,  &
171                     14.001 , 9.599 , 6.504 , 4.439 , 3.126 , 2.370 ,  &
172                     2.000/
173     
174    pressCLee = etaCLee*9.2e6
175   
176    DO j=jj_begin-1,jj_end+1
177       DO i=ii_begin-1,ii_end+1
178          ij=(j-1)*iim+i
179
180          DO l = 1, llm
181
182             pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff         
183             ! look for largest level such that pressCLee(level) > pplay(ij,l)) 
184             ! => pressClee(level+1) < pplay(ij,l) < pressClee(level)
185             level = 1
186             DO ll = 1, 30 ! 30 data levels
187                IF(pressCLee(ll) > pplay) THEN
188                   level = ll
189                END IF
190             END DO
191             
192             ! interpolate between level and level+1
193             ! interpolation is linear in log(pressure)
194             fact  = ( log10(pplay)-log10(pressCLee(level))) &
195                    /( log10(pressCLee(level+1))-log10(pressCLee(level)) )
196             ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact
197             zdt   = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact
198
199             IF(gradient) THEN
200                theta_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.)
201             ELSE
202                theta_eq(ij,l) = ztemp
203             END IF
204          END DO
205       END DO
206    END DO
207  END SUBROUTINE compute_temp_ref
208
209END MODULE etat0_venus_mod
Note: See TracBrowser for help on using the repository browser.