source: codes/icosagcm/trunk/src/etat0_dcmip5.f90 @ 186

Last change on this file since 186 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: 9.4 KB
Line 
1MODULE etat0_dcmip5_mod
2  USE icosa
3  PRIVATE
4  REAL(rstd),PARAMETER :: zt=15000
5  REAL(rstd),PARAMETER :: q0=0.021
6  REAL(rstd),PARAMETER :: qt=1e-11
7  REAL(rstd),PARAMETER :: T0=302.15
8  REAL(rstd),PARAMETER :: Ts=302.15
9  REAL(rstd),PARAMETER :: zq1=3000
10  REAL(rstd),PARAMETER :: zq2=8000
11  REAL(rstd),PARAMETER :: Gamma=0.007
12  REAL(rstd),PARAMETER :: pb=101500
13  REAL(rstd),PARAMETER :: Deltap=1115
14  REAL(rstd),PARAMETER :: rp=282000
15  REAL(rstd),PARAMETER :: zp=7000
16  REAL(rstd),PARAMETER :: epsilon=1e-25
17  REAL(rstd),PARAMETER :: Rd=287
18  REAL(rstd),SAVE :: lonc
19!$OMP THREADPRIVATE(lonc) 
20  REAL(rstd),SAVE :: latc
21!$OMP THREADPRIVATE(latc) 
22  TYPE(t_field),POINTER :: f_out_i(:)
23  REAL(rstd),POINTER,SAVE :: out_i(:,:)
24!$OMP THREADPRIVATE(out_i) 
25
26  PUBLIC  etat0
27CONTAINS
28 
29   
30   
31  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
32  USE icosa
33  IMPLICIT NONE
34    TYPE(t_field),POINTER :: f_ps(:)
35    TYPE(t_field),POINTER :: f_phis(:)
36    TYPE(t_field),POINTER :: f_theta_rhodz(:)
37    TYPE(t_field),POINTER :: f_u(:)
38    TYPE(t_field),POINTER :: f_q(:)
39 
40    REAL(rstd),POINTER :: ps(:)
41    REAL(rstd),POINTER :: phis(:)
42    REAL(rstd),POINTER :: theta_rhodz(:,:)
43    REAL(rstd),POINTER :: u(:,:)
44    REAL(rstd),POINTER :: q(:,:,:)
45    INTEGER :: ind
46 
47    CALL allocate_field(f_out_i,field_t,type_real,llm)
48 
49    DO ind=1,ndomain
50      IF (.NOT. assigned_domain(ind)) CYCLE
51      CALL swap_dimensions(ind)
52      CALL swap_geometry(ind)
53      ps=f_ps(ind)
54      phis=f_phis(ind)
55      theta_rhodz=f_theta_rhodz(ind)
56      u=f_u(ind)
57      q=f_q(ind)
58      out_i=f_out_i(ind)
59      CALL compute_etat0_dcmip5(ps, phis, theta_rhodz, u, q)
60    ENDDO
61    CALL writefield("out_i",f_out_i)
62    CALL deallocate_field(f_out_i)
63  END SUBROUTINE etat0
64 
65  SUBROUTINE compute_etat0_dcmip5(ps, phis, theta_rhodz, ue, q)
66  USE icosa
67  USE pression_mod
68  USE wind_mod
69  USE theta2theta_rhodz_mod
70  USE exner_mod
71  IMPLICIT NONE 
72  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
73  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
74  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
75  REAL(rstd),INTENT(OUT) :: ue(3*iim*jjm,llm)
76  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
77 
78  INTEGER :: i,j,l,ij
79  INTEGER :: testcase
80  REAL(rstd) :: Tv0, Tvt 
81  REAL(rstd) :: r 
82  REAL(rstd) :: z(iim*jjm,llm),zz 
83  REAL(rstd) :: p(iim*jjm,llm+1) 
84  REAL(rstd) :: Tave,Tp 
85  REAL(rstd) :: T(iim*jjm,llm) 
86  REAL(rstd) :: vt 
87  REAL(rstd) :: d1,d2,d 
88  REAL(rstd) :: ulon(3*iim*jjm,llm) 
89  REAL(rstd) :: ulat(3*iim*jjm,llm)
90  REAL(rstd) :: lon,lat
91  REAL(rstd) :: pks(iim*jjm) 
92  REAL(rstd) :: pk(iim*jjm,llm) 
93
94
95    latc=Pi/18
96    lonc=Pi 
97
98    Tv0=T0/(1+0.608*q0)
99    Tvt=Tv0-Gamma*zt
100 
101    testcase=1
102    CALL getin("dcmip5_testcase",testcase)
103   
104    phis(:)=0.
105   
106    DO j=jj_begin,jj_end
107      DO i=ii_begin,ii_end
108        ij=(j-1)*iim+i
109        CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
110        r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
111        ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5)
112      ENDDO
113    ENDDO
114   
115    CALL compute_pression(ps,p,0)
116
117    DO l=1,llm
118      DO j=jj_begin,jj_end
119        DO i=ii_begin,ii_end
120          ij=(j-1)*iim+i
121          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
122          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
123          CALL compute_z(0.5*(p(ij,l)+p(ij,l+1)),ps(ij),r,z(ij,l))
124        ENDDO
125      ENDDO
126    ENDDO
127    out_i=z
128   
129    DO l=1,llm
130      DO j=jj_begin-1,jj_end+1
131        DO i=ii_begin-1,ii_end+1
132          ij=(j-1)*iim+i
133          IF (z(ij,l)<=zt) THEN
134            q(ij,l,1)=q0*exp(-z(ij,l)/zq1)*exp(-(z(ij,l)/zq2)**2)
135          ELSE
136            q(ij,l,1)=qt
137          ENDIF
138        ENDDO
139      ENDDO
140    ENDDO
141
142   
143     DO l=1,llm
144      DO j=jj_begin-1,jj_end+1
145        DO i=ii_begin-1,ii_end+1
146          ij=(j-1)*iim+i
147          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
148          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
149          IF (z(ij,l)<=zt) THEN
150            Tave=Tv0-Gamma*z(ij,l)
151            Tp=(Tv0-Gamma*z(ij,l)) * ( ( 1 + 2*Rd*(Tv0-Gamma*z(ij,l))*z(ij,l) /                           &
152                                            (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(z(ij,l)/zp)**2))))**(-1) - 1) 
153          ELSE
154            Tave=Tvt
155            Tp=0
156          ENDIF
157          T(ij,l)=Tave+Tp 
158          out_i(ij,l)=Tave+Tp
159        ENDDO
160      ENDDO
161    ENDDO
162
163    DO l=1,llm
164      DO j=jj_begin-1,jj_end+1
165        DO i=ii_begin-1,ii_end+1
166          ij=(j-1)*iim+i
167
168          CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
169          zz=0.5*(z(ij,l)+z(ij+t_right,l))
170          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
171          IF (zz<=zt) THEN
172            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
173                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
174                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
175          ELSE
176            vt = 0
177          ENDIF
178          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
179          d2=cos(latc)*sin(lon-lonc)
180          d=MAX(1e-25,sqrt(d1**2+d2**2))
181          ulon(ij+u_right,l)=vt*d1/d
182          ulat(ij+u_right,l)=vt*d2/d
183
184
185         
186          CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
187          zz=0.5*(z(ij,l)+z(ij+t_lup,l))
188          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
189          IF (zz<=zt) THEN
190            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
191                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
192                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
193
194         ELSE
195            vt = 0
196          ENDIF
197          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
198          d2=cos(latc)*sin(lon-lonc)
199          d=MAX(1e-25,sqrt(d1**2+d2**2))
200          ulon(ij+u_lup,l)=vt*d1/d
201          ulat(ij+u_lup,l)=vt*d2/d
202         
203         
204         
205         
206          CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
207          zz=0.5*(z(ij,l)+z(ij+t_ldown,l))
208          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
209          IF (zz<=zt) THEN
210            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
211                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
212                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
213
214          ELSE
215            vt = 0
216          ENDIF
217          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
218          d2=cos(latc)*sin(lon-lonc)
219          d=MAX(1e-25,sqrt(d1**2+d2**2))
220          ulon(ij+u_ldown,l)=vt*d1/d
221          ulat(ij+u_ldown,l)=vt*d2/d
222         
223
224
225          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
226          zz=z(ij,l)
227          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
228          IF (zz<=zt) THEN
229            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
230                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
231                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
232          ELSE
233            vt = 0
234          ENDIF
235          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
236          d2=cos(latc)*sin(lon-lonc)
237          d=MAX(1e-25,sqrt(d1**2+d2**2))
238!          ulon(ij+u_ldown,l)=vt*d1/d
239!          ulat(ij+u_ldown,l)=vt*d2/d
240!          out_i(ij,l)=sqrt((vt*d1/d)**2+(vt*d2/d)**2)
241        ENDDO
242      ENDDO
243    ENDDO
244   
245    CALL compute_wind_perp_from_lonlat_compound(ulon, ulat, ue)
246    CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
247
248  END SUBROUTINE compute_etat0_dcmip5
249 
250  SUBROUTINE compute_z(pmodel,ps,r,z)
251  USE icosa
252  IMPLICIT NONE
253     REAL(rstd),PARAMETER   :: epsilon0=2e-13
254     REAL(rstd),INTENT(IN)  :: pmodel
255     REAL(rstd),INTENT(IN)  :: ps
256     REAL(rstd),INTENT(IN)  :: r
257     REAL(rstd),INTENT(OUT) :: z
258     
259     REAL(rstd) :: Tv0,Tvt,pt
260     REAL(rstd) :: pave, pp, dpdz
261     REAL(rstd) :: znp1
262     REAL(rstd) :: epsilon
263     REAL(rstd) :: p
264     INTEGER  :: n
265         
266     Tv0=T0/(1+0.608*q0)
267     Tvt=Tv0-Gamma*zt
268     pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma))
269
270     IF (pmodel>pt) THEN
271       z=Tv0/Gamma*(1-(pmodel/ps)**(Rd*Gamma/g))
272     ELSE
273       z=zt+Rd*Tvt/g*log(Pt/pmodel)
274     ENDIF
275     
276!     PRINT*,pmodel,pt,z,zt
277     IF (z>zt .OR. r>1000000) return
278     
279     epsilon=1
280     n=0
281     
282     DO WHILE(epsilon>epsilon0 .AND. n<20)
283 
284       IF (z<zt) THEN
285         pave=pb*((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))
286         pp= -Deltap * exp(-(r/rp)**1.5-(z/zp)**2) *  ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma))
287       ELSE
288         pave=pt*exp(g*(zt-z)/(Rd*Tvt))
289         pp=0
290       ENDIF
291       p=pave+pp
292     
293       dpdz =  2*Deltap*z/zp**2 * exp(-(r/rp)**1.5) * exp(-(z/zp)**2) * ((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))  &
294             - g/(Rd*Tv0) * (preff - Deltap*exp(-(r/rp)**1.5) * exp(-(z/zp)**2)) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)-1)
295           
296       znp1 = z - ( pmodel - p ) / (-dpdz)
297       
298       epsilon=ABS( (znp1-z)/znp1 )
299!       PRINT *,"epsilon",epsilon,r,z,znp1
300!       PRINT *,"----",n,pmodel,ps,pave,dpdz
301       z=znp1
302       n=n+1
303     ENDDO 
304   
305   END SUBROUTINE compute_z
306
307END MODULE etat0_dcmip5_mod
Note: See TracBrowser for help on using the repository browser.