source: dynamico_lmdz/aquaplanet/ICOSAGCM/src/etat0_venus.f90 @ 3891

Last change on this file since 3891 was 3810, checked in by ymipsl, 10 years ago

Add DYNAMICO in aquaplanet configuration
YM
  • Property svn:executable set to *
File size: 7.0 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    REAL(rstd) :: lat(iim*jjm)        ! latitude                   
106    REAL(rstd) :: pplay(iim*jjm, llm) ! pressure at full layers
107    INTEGER :: ind
108
109    CALL allocate_field(f_temp,field_t,type_real,llm)
110    DO ind=1,ndomain
111       IF (.NOT. assigned_domain(ind)) CYCLE
112       CALL swap_dimensions(ind)
113       CALL swap_geometry(ind)
114       ps=f_ps(ind)
115       ps(:)=preff
116       phis=f_phis(ind)
117       phis(:)=0.
118       u=f_u(ind)
119       u(:,:)=0
120       q=f_q(ind)
121       q(:,:,:)=1e2
122       temp=f_temp(ind)
123       CALL compute_temp_ref(temp, .FALSE.) ! Without meridional gradient
124    ENDDO
125    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
126    CALL deallocate_field(f_temp)
127
128  END SUBROUTINE etat0
129
130!------------------------- Compute reference temperature field ------------------------
131
132  SUBROUTINE compute_temp_ref(theta_eq, gradient)
133    USE icosa
134    USE disvert_mod
135    USE omp_para
136    USE math_const
137    IMPLICIT NONE
138    REAL(rstd), INTENT(OUT) :: theta_eq(iim*jjm,llm)
139    LOGICAL, INTENT(IN) :: gradient
140    REAL(rstd)  :: clat(iim*jjm)        ! latitude
141    integer :: level
142    REAL ::  pressCLee(31), tempCLee(31), dt_epCLee(31), etaCLee(31)
143   
144    real(rstd) ::  lon,lat, pplay, ztemp,zdt,fact
145    logical, save ::  firstcall
146    integer :: i,j,ij, l,ll
147   
148    data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, &
149                      4.273e-1, 3.373e-1, 2.610e-1,1.979e-1,1.472e-1,  &
150                      1.074e-1, 7.672e-2, 5.361e-2,3.657e-2,2.430e-2,  &
151                      1.569e-2, 9.814e-3, 5.929e-3,3.454e-3,1.934e-3,  &
152                      1.043e-3, 5.400e-4, 2.710e-4,1.324e-4,6.355e-5,  &
153                      3.070e-5, 1.525e-5, 7.950e-6,4.500e-6,2.925e-6,  &
154                      2.265e-6/
155
156    DO j=jj_begin-1,jj_end+1
157       DO i=ii_begin-1,ii_end+1
158          ij=(j-1)*iim+i
159          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
160          clat(ij)=cos(lat)
161       ENDDO
162    ENDDO
163
164    data   tempCLee/ 728.187, 715.129, 697.876, 677.284, 654.078, 628.885, &
165                     602.225, 574.542, 546.104, 517.339, 488.560, 459.932, &
166                     431.741, 404.202, 377.555, 352.042, 327.887, 305.313, &
167                     284.556, 265.697, 248.844, 233.771, 220.368, 208.247, &
168                     197.127, 187.104, 178.489, 171.800, 167.598, 165.899, &
169                     165.676/
170    data   dt_epCLee/6.101 , 6.136 , 6.176 , 6.410 , 6.634 , 6.678 , &
171                     6.719 , 6.762 , 7.167 , 7.524 , 9.840 ,14.948 ,  &
172                     21.370 ,28.746 ,36.373 ,43.315 ,48.534 ,51.175 ,  &
173                     50.757 ,47.342 ,41.536 ,34.295 ,26.758 ,19.807 ,  &
174                     14.001 , 9.599 , 6.504 , 4.439 , 3.126 , 2.370 ,  &
175                     2.000/
176     
177    pressCLee = etaCLee*9.2e6
178   
179    DO j=jj_begin-1,jj_end+1
180       DO i=ii_begin-1,ii_end+1
181          ij=(j-1)*iim+i
182
183          DO l = 1, llm
184
185             pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff         
186             ! look for largest level such that pressCLee(level) > pplay(ij,l)) 
187             ! => pressClee(level+1) < pplay(ij,l) < pressClee(level)
188             level = 1
189             DO ll = 1, 30 ! 30 data levels
190                IF(pressCLee(ll) > pplay) THEN
191                   level = ll
192                END IF
193             END DO
194             
195             ! interpolate between level and level+1
196             ! interpolation is linear in log(pressure)
197             fact  = ( log10(pplay)-log10(pressCLee(level))) &
198                    /( log10(pressCLee(level+1))-log10(pressCLee(level)) )
199             ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact
200             zdt   = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact
201
202             IF(gradient) THEN
203                theta_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.)
204             ELSE
205                theta_eq(ij,l) = ztemp
206             END IF
207          END DO
208       END DO
209    END DO
210  END SUBROUTINE compute_temp_ref
211
212END MODULE etat0_venus_mod
Note: See TracBrowser for help on using the repository browser.