source: codes/icosagcm/trunk/src/etat0_jablonowsky06.f90 @ 17

Last change on this file since 17 was 17, checked in by ymipsl, 12 years ago

Merge advection scheme from sarvesh in standard version

YM
File size: 7.6 KB
Line 
1MODULE etat0_jablonowsky06_mod
2  USE genmod
3  PRIVATE
4  REAL(rstd),PARAMETER :: eta0=0.252
5  REAL(rstd),PARAMETER :: etat=0.2
6  REAL(rstd),PARAMETER :: ps0=1e5
7  REAL(rstd),PARAMETER :: u0=35
8!  REAL(rstd),PARAMETER :: u0=0
9  REAL(rstd),PARAMETER :: T0=288
10  REAL(rstd),PARAMETER :: DeltaT=4.8e5
11  REAL(rstd),PARAMETER :: Rd=287
12  REAL(rstd),PARAMETER :: Gamma=0.005
13  REAL(rstd),PARAMETER :: up0=1
14  PUBLIC  test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06
15CONTAINS
16 
17  SUBROUTINE test_etat0_jablonowsky06
18  USE field_mod
19  USE domain_mod
20  USE dimensions
21  USE grid_param
22  USE geometry
23  USE write_field
24  USE kinetic_mod
25  USE pression_mod
26  USE exner_mod
27  USE geopotential_mod
28  USE vorticity_mod
29  IMPLICIT NONE
30    TYPE(t_field),POINTER :: f_ps(:)
31    TYPE(t_field),POINTER :: f_phis(:)
32    TYPE(t_field),POINTER :: f_theta_rhodz(:)
33    TYPE(t_field),POINTER :: f_u(:)
34    TYPE(t_field),POINTER :: f_Ki(:)
35    TYPE(t_field),POINTER :: f_temp(:)
36    TYPE(t_field),POINTER :: f_p(:)
37    TYPE(t_field),POINTER :: f_pks(:)
38    TYPE(t_field),POINTER :: f_pk(:)
39    TYPE(t_field),POINTER :: f_phi(:)
40    TYPE(t_field),POINTER :: f_vort(:)
41    TYPE(t_field),POINTER :: f_q(:)
42 
43    REAL(rstd),POINTER :: Ki(:,:)
44    REAL(rstd),POINTER :: temp(:)
45    INTEGER :: ind
46   
47   
48    CALL allocate_field(f_ps,field_t,type_real)
49    CALL allocate_field(f_phis,field_t,type_real)
50    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
51    CALL allocate_field(f_u,field_u,type_real,llm)
52    CALL allocate_field(f_p,field_t,type_real,llm+1)
53    CALL allocate_field(f_Ki,field_t,type_real,llm)
54    CALL allocate_field(f_pks,field_t,type_real)
55    CALL allocate_field(f_pk,field_t,type_real,llm)
56    CALL allocate_field(f_phi,field_t,type_real,llm)
57    CALL allocate_field(f_temp,field_t,type_real)
58    CALL allocate_field(f_vort,field_z,type_real,llm)
59   
60    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
61
62    CALL kinetic(f_u,f_Ki)
63    CALL vorticity(f_u,f_vort)
64
65    CALL pression(f_ps,f_p)     
66    CALL exner(f_ps,f_p,f_pks,f_pk) 
67    CALL geopotential(f_phis,f_pks,f_pk,f_theta_rhodz,f_phi)
68
69    CALL writefield('ps',f_ps)
70    CALL writefield('phis',f_phis)
71    CALL writefield('theta',f_theta_rhodz)
72    CALL writefield('f_phi',f_phi)
73
74    CALL writefield('Ki',f_Ki)
75    CALL writefield('vort',f_vort)
76   
77  END SUBROUTINE test_etat0_jablonowsky06
78   
79   
80  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
81  USE field_mod
82  USE domain_mod
83  USE domain_mod
84  USE dimensions
85  USE grid_param
86  USE geometry
87  IMPLICIT NONE
88    TYPE(t_field),POINTER :: f_ps(:)
89    TYPE(t_field),POINTER :: f_phis(:)
90    TYPE(t_field),POINTER :: f_theta_rhodz(:)
91    TYPE(t_field),POINTER :: f_u(:)
92    TYPE(t_field),POINTER :: f_q(:)
93 
94    REAL(rstd),POINTER :: ps(:)
95    REAL(rstd),POINTER :: phis(:)
96    REAL(rstd),POINTER :: theta_rhodz(:,:)
97    REAL(rstd),POINTER :: u(:,:)
98    REAL(rstd),POINTER :: q(:,:,:)
99    INTEGER :: ind
100   
101    DO ind=1,ndomain
102      CALL swap_dimensions(ind)
103      CALL swap_geometry(ind)
104      ps=f_ps(ind)
105      phis=f_phis(ind)
106      theta_rhodz=f_theta_rhodz(ind)
107      u=f_u(ind)
108      q=f_q(ind)
109      q=0
110      CALL compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u)
111    ENDDO
112
113  END SUBROUTINE etat0
114 
115  SUBROUTINE compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u)
116  USE domain_mod
117  USE dimensions
118  USE grid_param
119  USE geometry
120  USE metric
121  USE disvert_mod
122  USE spherical_geom_mod
123  USE vector
124  USE pression_mod
125  USE exner_mod
126  USE geopotential_mod
127  USE theta2theta_rhodz_mod
128  IMPLICIT NONE 
129  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
130  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
131  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
132  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
133 
134  INTEGER :: i,j,l,ij
135  REAL(rstd) :: theta(iim*jjm,llm)
136  REAL(rstd) :: eta(llm)
137  REAL(rstd) :: etav(llm)
138  REAL(rstd) :: etas, etavs
139  REAL(rstd) :: lon,lat
140  REAL(rstd) :: ulon(3)
141  REAL(rstd) :: ep(3), norm_ep
142  REAL(rstd) :: Tave, T
143  REAL(rstd) :: phis_ave
144  REAL(rstd) :: V0(3)
145  REAL(rstd) :: r2
146  REAL(rstd) :: utot
147  REAL(rstd) :: lonx,latx
148 
149    DO l=1,llm
150      eta(l)= 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
151      etav(l)=(eta(l)-eta0)*Pi/2
152    ENDDO
153    etas=ap(1)+bp(1)
154    etavs=(etas-eta0)*Pi/2
155
156    DO j=jj_begin,jj_end
157      DO i=ii_begin,ii_end
158        ij=(j-1)*iim+i
159        ps(ij)=ps0
160      ENDDO
161    ENDDO
162   
163   
164    CALL lonlat2xyz(Pi/9,2*Pi/9,V0)
165
166    u(:,:)=1e10     
167    DO l=1,llm
168      DO j=jj_begin,jj_end
169        DO i=ii_begin,ii_end
170          ij=(j-1)*iim+i
171         
172          CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat)
173          CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep)
174          r2=(asin(sqrt(sum(ep*ep))))**2
175          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
176
177          ulon(1) = -sin(lon) * utot
178          ulon(2) = cos(lon)  * utot
179          ulon(3) = 0         * utot
180         
181                   
182          CALL cross_product2(xyz_v(ij+z_rdown,:)/radius,xyz_v(ij+z_rup,:)/radius,ep)
183          norm_ep=sqrt(sum(ep(:)**2))
184          IF (norm_ep>1e-30) THEN
185            ep=-ep/norm_ep*ne(ij,right)
186            u(ij+u_right,l)=sum(ep(:)*ulon(:))
187          ENDIF
188
189
190         
191          CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat)
192          CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep)
193          r2=(asin(sqrt(sum(ep*ep))))**2
194          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
195          ulon(1) = -sin(lon) * utot
196          ulon(2) = cos(lon)  * utot
197          ulon(3) = 0         * utot
198         
199                   
200          CALL cross_product2(xyz_v(ij+z_up,:)/radius,xyz_v(ij+z_lup,:)/radius,ep)
201          norm_ep=sqrt(sum(ep(:)**2))
202          IF (norm_ep>1e-30) THEN
203            ep=-ep/norm_ep*ne(ij,lup)
204            u(ij+u_lup,l)=sum(ep(:)*ulon(:))
205          ENDIF
206
207
208
209
210
211                   
212          CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat)
213          CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep)
214          r2=(asin(sqrt(sum(ep*ep))))**2
215          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
216          ulon(1) = -sin(lon) * utot
217          ulon(2) = cos(lon)  * utot
218          ulon(3) = 0         * utot
219         
220                   
221          CALL cross_product2(xyz_v(ij+z_ldown,:)/radius,xyz_v(ij+z_down,:)/radius,ep)
222          norm_ep=sqrt(sum(ep(:)**2))
223          IF (norm_ep>1e-30) THEN
224            ep=-ep/norm_ep*ne(ij,ldown)
225            u(ij+u_ldown,l)=sum(ep(:)*ulon(:))
226          ENDIF         
227
228       ENDDO
229     ENDDO
230    ENDDO 
231     
232     
233     DO l=1,llm
234       Tave=T0*eta(l)**(Rd*Gamma/g)
235       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
236       DO j=jj_begin,jj_end
237         DO i=ii_begin,ii_end
238           ij=(j-1)*iim+i
239           CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
240             
241            T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5              &
242                             * ( (-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
243                                + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega)
244           
245            theta(ij,l)=T*eta(l)**(-kappa)
246
247          ENDDO
248       ENDDO
249     ENDDO
250     
251     
252     phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
253     DO j=jj_begin,jj_end
254       DO i=ii_begin,ii_end
255         ij=(j-1)*iim+i
256         CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
257         phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
258                                           +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega )
259       ENDDO
260     ENDDO
261
262    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
263
264  END SUBROUTINE compute_etat0_jablonowsky06
265 
266END MODULE etat0_jablonowsky06_mod
Note: See TracBrowser for help on using the repository browser.