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

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

dynamico tree creation

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 
143    DO l=1,llm
144      eta(l)= 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
145      etav(l)=(eta(l)-eta0)*Pi/2
146    ENDDO
147    etas=ap(1)+bp(1)
148    etavs=(etas-eta0)*Pi/2
149
150    DO j=jj_begin,jj_end
151      DO i=ii_begin,ii_end
152        ij=(j-1)*iim+i
153        ps(ij)=ps0
154      ENDDO
155    ENDDO
156   
157   
158    CALL lonlat2xyz(Pi/9,2*Pi/9,V0)
159    u(:,:)=1e10     
160    DO l=1,llm
161      DO j=jj_begin,jj_end
162        DO i=ii_begin,ii_end
163          ij=(j-1)*iim+i
164         
165          CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat)
166          CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep)
167          r2=(asin(sum(ep*ep)))**2
168          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
169
170          ulon(1) = -sin(lon) * utot
171          ulon(2) = cos(lon)  * utot
172          ulon(3) = 0         * utot
173         
174                   
175          CALL cross_product2(xyz_v(ij+z_rdown,:)/radius,xyz_v(ij+z_rup,:)/radius,ep)
176          norm_ep=sqrt(sum(ep(:)**2))
177          IF (norm_ep>1e-30) THEN
178            ep=-ep/norm_ep*ne(ij,right)
179            u(ij+u_right,l)=sum(ep(:)*ulon(:))
180          ENDIF
181
182
183         
184          CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat)
185          CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep)
186          r2=(asin(sum(ep*ep)))**2
187          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
188          ulon(1) = -sin(lon) * utot
189          ulon(2) = cos(lon)  * utot
190          ulon(3) = 0         * utot
191         
192                   
193          CALL cross_product2(xyz_v(ij+z_up,:)/radius,xyz_v(ij+z_lup,:)/radius,ep)
194          norm_ep=sqrt(sum(ep(:)**2))
195          IF (norm_ep>1e-30) THEN
196            ep=-ep/norm_ep*ne(ij,lup)
197            u(ij+u_lup,l)=sum(ep(:)*ulon(:))
198          ENDIF
199
200
201
202
203
204                   
205          CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat)
206          CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep)
207          r2=(asin(sum(ep*ep)))**2
208          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
209          ulon(1) = -sin(lon) * utot
210          ulon(2) = cos(lon)  * utot
211          ulon(3) = 0         * utot
212         
213                   
214          CALL cross_product2(xyz_v(ij+z_ldown,:)/radius,xyz_v(ij+z_down,:)/radius,ep)
215          norm_ep=sqrt(sum(ep(:)**2))
216          IF (norm_ep>1e-30) THEN
217            ep=-ep/norm_ep*ne(ij,ldown)
218            u(ij+u_ldown,l)=sum(ep(:)*ulon(:))
219          ENDIF         
220
221       ENDDO
222     ENDDO
223    ENDDO 
224     
225     
226     DO l=1,llm
227       Tave=T0*eta(l)**(Rd*Gamma/g)
228       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
229       DO j=jj_begin,jj_end
230         DO i=ii_begin,ii_end
231           ij=(j-1)*iim+i
232           CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
233             
234            T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5              &
235                             * ( (-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
236                                + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega)
237           
238            theta(ij,l)=T*eta(l)**(-kappa)
239
240          ENDDO
241       ENDDO
242     ENDDO
243     
244     
245     phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
246     DO j=jj_begin,jj_end
247       DO i=ii_begin,ii_end
248         ij=(j-1)*iim+i
249         CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
250         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  &
251                                           +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega )
252       ENDDO
253     ENDDO
254
255    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
256
257  END SUBROUTINE compute_etat0_jablonowsky06
258 
259END MODULE etat0_jablonowsky06_mod
Note: See TracBrowser for help on using the repository browser.