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

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

Fixed DCMIP5 physics/etat0
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.