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