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

Last change on this file since 191 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM
File size: 7.4 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,SAVE :: f_ps(:)
26    TYPE(t_field),POINTER,SAVE :: f_phis(:)
27    TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:)
28    TYPE(t_field),POINTER,SAVE :: f_u(:)
29    TYPE(t_field),POINTER,SAVE :: f_Ki(:)
30    TYPE(t_field),POINTER,SAVE :: f_temp(:)
31    TYPE(t_field),POINTER,SAVE :: f_p(:)
32    TYPE(t_field),POINTER,SAVE :: f_pks(:)
33    TYPE(t_field),POINTER,SAVE :: f_pk(:)
34    TYPE(t_field),POINTER,SAVE :: f_phi(:)
35    TYPE(t_field),POINTER,SAVE :: f_vort(:)
36    TYPE(t_field),POINTER,SAVE :: 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      IF (.NOT. assigned_domain(ind)) CYCLE
93      CALL swap_dimensions(ind)
94      CALL swap_geometry(ind)
95      ps=f_ps(ind)
96      phis=f_phis(ind)
97      theta_rhodz=f_theta_rhodz(ind)
98      u=f_u(ind)
99      q=f_q(ind)
100      q=0
101      CALL compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u)
102    ENDDO
103
104  END SUBROUTINE etat0
105 
106  SUBROUTINE compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u)
107  USE icosa
108  USE disvert_mod
109  USE pression_mod
110  USE exner_mod
111  USE geopotential_mod
112  USE theta2theta_rhodz_mod
113  IMPLICIT NONE 
114  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
115  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
116  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
117  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
118 
119  INTEGER :: i,j,l,ij
120  REAL(rstd) :: theta(iim*jjm,llm)
121  REAL(rstd) :: eta(llm)
122  REAL(rstd) :: etav(llm)
123  REAL(rstd) :: etas, etavs
124  REAL(rstd) :: lon,lat
125  REAL(rstd) :: ulon(3)
126  REAL(rstd) :: ep(3), norm_ep
127  REAL(rstd) :: Tave, T
128  REAL(rstd) :: phis_ave
129  REAL(rstd) :: V0(3)
130  REAL(rstd) :: r2
131  REAL(rstd) :: utot
132  REAL(rstd) :: lonx,latx
133 
134    DO l=1,llm
135      eta(l)= 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
136      etav(l)=(eta(l)-eta0)*Pi/2
137    ENDDO
138    etas=ap(1)+bp(1)
139    etavs=(etas-eta0)*Pi/2
140
141    DO j=jj_begin,jj_end
142      DO i=ii_begin,ii_end
143        ij=(j-1)*iim+i
144        ps(ij)=ps0
145      ENDDO
146    ENDDO
147   
148   
149    CALL lonlat2xyz(Pi/9,2*Pi/9,V0)
150
151    u(:,:)=1e10     
152    DO l=1,llm
153      DO j=jj_begin,jj_end
154        DO i=ii_begin,ii_end
155          ij=(j-1)*iim+i
156         
157          CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat)
158          CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep)
159          r2=(asin(sqrt(sum(ep*ep))))**2
160          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
161
162          ulon(1) = -sin(lon) * utot
163          ulon(2) = cos(lon)  * utot
164          ulon(3) = 0         * utot
165         
166                   
167          CALL cross_product2(xyz_v(ij+z_rdown,:)/radius,xyz_v(ij+z_rup,:)/radius,ep)
168          norm_ep=sqrt(sum(ep(:)**2))
169          IF (norm_ep>1e-30) THEN
170            ep=-ep/norm_ep*ne(ij,right)
171            u(ij+u_right,l)=sum(ep(:)*ulon(:))
172          ENDIF
173
174
175         
176          CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat)
177          CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep)
178          r2=(asin(sqrt(sum(ep*ep))))**2
179          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
180          ulon(1) = -sin(lon) * utot
181          ulon(2) = cos(lon)  * utot
182          ulon(3) = 0         * utot
183         
184                   
185          CALL cross_product2(xyz_v(ij+z_up,:)/radius,xyz_v(ij+z_lup,:)/radius,ep)
186          norm_ep=sqrt(sum(ep(:)**2))
187          IF (norm_ep>1e-30) THEN
188            ep=-ep/norm_ep*ne(ij,lup)
189            u(ij+u_lup,l)=sum(ep(:)*ulon(:))
190          ENDIF
191
192
193
194
195
196                   
197          CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat)
198          CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep)
199          r2=(asin(sqrt(sum(ep*ep))))**2
200          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01)
201          ulon(1) = -sin(lon) * utot
202          ulon(2) = cos(lon)  * utot
203          ulon(3) = 0         * utot
204         
205                   
206          CALL cross_product2(xyz_v(ij+z_ldown,:)/radius,xyz_v(ij+z_down,:)/radius,ep)
207          norm_ep=sqrt(sum(ep(:)**2))
208          IF (norm_ep>1e-30) THEN
209            ep=-ep/norm_ep*ne(ij,ldown)
210            u(ij+u_ldown,l)=sum(ep(:)*ulon(:))
211          ENDIF         
212
213       ENDDO
214     ENDDO
215    ENDDO 
216     
217     
218     DO l=1,llm
219       Tave=T0*eta(l)**(Rd*Gamma/g)
220       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
221       DO j=jj_begin,jj_end
222         DO i=ii_begin,ii_end
223           ij=(j-1)*iim+i
224           CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
225             
226            T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5              &
227                             * ( (-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
228                                + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega)
229           
230            theta(ij,l)=T*eta(l)**(-kappa)
231
232          ENDDO
233       ENDDO
234     ENDDO
235     
236     
237     phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
238     DO j=jj_begin,jj_end
239       DO i=ii_begin,ii_end
240         ij=(j-1)*iim+i
241         CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
242         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  &
243                                           +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega )
244       ENDDO
245     ENDDO
246
247    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
248
249  END SUBROUTINE compute_etat0_jablonowsky06
250 
251END MODULE etat0_jablonowsky06_mod
Note: See TracBrowser for help on using the repository browser.