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

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

Simplify the management of the module.

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