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

Last change on this file since 201 was 201, checked in by dubos, 10 years ago

Simplified etat0 : isothermal, jablonowsky06
File size: 9.1 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 :: T0=288
9  REAL(rstd),PARAMETER :: DeltaT=4.8e5
10  REAL(rstd),PARAMETER :: Rd=287
11  REAL(rstd),PARAMETER :: Gamma=0.005
12  REAL(rstd),PARAMETER :: up0=1
13  PUBLIC  test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06, compute_etat0_new
14CONTAINS
15 
16  SUBROUTINE test_etat0_jablonowsky06
17  USE icosa
18  USE kinetic_mod
19  USE pression_mod
20  USE exner_mod
21  USE geopotential_mod
22  USE vorticity_mod
23  IMPLICIT NONE
24    TYPE(t_field),POINTER,SAVE :: f_ps(:)
25    TYPE(t_field),POINTER,SAVE :: f_phis(:)
26    TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:)
27    TYPE(t_field),POINTER,SAVE :: f_u(:)
28    TYPE(t_field),POINTER,SAVE :: f_Ki(:)
29    TYPE(t_field),POINTER,SAVE :: f_temp(:)
30    TYPE(t_field),POINTER,SAVE :: f_p(:)
31    TYPE(t_field),POINTER,SAVE :: f_pks(:)
32    TYPE(t_field),POINTER,SAVE :: f_pk(:)
33    TYPE(t_field),POINTER,SAVE :: f_phi(:)
34    TYPE(t_field),POINTER,SAVE :: f_vort(:)
35    TYPE(t_field),POINTER,SAVE :: f_q(:)
36 
37    REAL(rstd),POINTER :: Ki(:,:)
38    REAL(rstd),POINTER :: temp(:)
39    INTEGER :: ind
40   
41   
42    CALL allocate_field(f_ps,field_t,type_real)
43    CALL allocate_field(f_phis,field_t,type_real)
44    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
45    CALL allocate_field(f_u,field_u,type_real,llm)
46    CALL allocate_field(f_p,field_t,type_real,llm+1)
47    CALL allocate_field(f_Ki,field_t,type_real,llm)
48    CALL allocate_field(f_pks,field_t,type_real)
49    CALL allocate_field(f_pk,field_t,type_real,llm)
50    CALL allocate_field(f_phi,field_t,type_real,llm)
51    CALL allocate_field(f_temp,field_t,type_real)
52    CALL allocate_field(f_vort,field_z,type_real,llm)
53   
54    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
55
56    CALL kinetic(f_u,f_Ki)
57    CALL vorticity(f_u,f_vort)
58
59    CALL pression(f_ps,f_p)     
60    CALL exner(f_ps,f_p,f_pks,f_pk) 
61    CALL geopotential(f_phis,f_pks,f_pk,f_theta_rhodz,f_phi)
62
63    CALL writefield('ps',f_ps)
64    CALL writefield('phis',f_phis)
65    CALL writefield('theta',f_theta_rhodz)
66    CALL writefield('f_phi',f_phi)
67
68    CALL writefield('Ki',f_Ki)
69    CALL writefield('vort',f_vort)
70   
71  END SUBROUTINE test_etat0_jablonowsky06
72   
73   
74  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
75  USE icosa
76  IMPLICIT NONE
77    TYPE(t_field),POINTER :: f_ps(:)
78    TYPE(t_field),POINTER :: f_phis(:)
79    TYPE(t_field),POINTER :: f_theta_rhodz(:)
80    TYPE(t_field),POINTER :: f_u(:)
81    TYPE(t_field),POINTER :: f_q(:)
82 
83    REAL(rstd),POINTER :: ps(:)
84    REAL(rstd),POINTER :: phis(:)
85    REAL(rstd),POINTER :: theta_rhodz(:,:)
86    REAL(rstd),POINTER :: u(:,:)
87    REAL(rstd),POINTER :: q(:,:,:)
88    INTEGER :: ind
89   
90    DO ind=1,ndomain
91      IF (.NOT. assigned_domain(ind)) CYCLE
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
250  SUBROUTINE compute_etat0_new(ngrid,lat,lon, phis, ps, temp, ulon, ulat)
251    USE disvert_mod
252    IMPLICIT NONE 
253    INTEGER, INTENT(IN) :: ngrid
254    REAL(rstd),INTENT(IN) :: lat(ngrid)
255    REAL(rstd),INTENT(IN) :: lon(ngrid)
256    REAL(rstd),INTENT(OUT) :: phis(ngrid)
257    REAL(rstd),INTENT(OUT) :: ps(ngrid)
258    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
259    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
260    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
261   
262    INTEGER :: l,ij
263    REAL(rstd) :: eta(llm)
264    REAL(rstd) :: etav(llm)
265    REAL(rstd) :: etas, etavs, Tave, phis_ave, r2
266   
267    DO l=1,llm
268       eta(l)  = 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
269       etav(l) = (eta(l)-eta0)*Pi/2
270    ENDDO
271    etas  = ap(1)+bp(1)
272    etavs = (etas-eta0)*Pi/2
273       
274    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
275    DO ij=1,ngrid
276       ps(ij)=ps0
277       phis(ij) = phis_ave + u0*cos(etavs)**1.5*( (-2*sin(lat(ij))**6 * (cos(lat(ij))**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
278            +(8./5*cos(lat(ij))**3 * (sin(lat(ij))**2 + 2./3) - Pi/4)*radius*Omega )
279    ENDDO
280
281    DO l=1,llm
282       Tave=T0*eta(l)**(Rd*Gamma/g)
283       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
284       DO ij=1,ngrid
285          r2 = arc(Pi/9,2*Pi/9, lon(ij),lat(ij))**2
286          temp(ij,l) = Tave + 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5  &
287               * ( (-2*sin(lat(ij))**6*(cos(lat(ij))**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
288               + (8./5*cos(lat(ij))**3*(sin(lat(ij))**2+2./3)-Pi/4)*radius*Omega)
289          ulon(ij,l) = u0*cos(etav(l))**1.5*sin(2*lat(ij))**2 + up0*exp(-r2/0.01)
290          ulat(ij,l) = 0
291       ENDDO
292    ENDDO
293   
294   
295  END SUBROUTINE compute_etat0_new
296 
297END MODULE etat0_jablonowsky06_mod
Note: See TracBrowser for help on using the repository browser.