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

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

Update on 3D dynamic

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