1 | MODULE etat0_dcmip4_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 | REAL(rstd) :: latw=2*Pi/9 |
---|
14 | REAL(rstd) :: pw=34000 |
---|
15 | REAL(rstd) :: q0=0.021 |
---|
16 | REAL(rstd) :: lonc |
---|
17 | REAL(rstd) :: latc |
---|
18 | PUBLIC etat0 |
---|
19 | CONTAINS |
---|
20 | |
---|
21 | |
---|
22 | |
---|
23 | SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) |
---|
24 | USE icosa |
---|
25 | IMPLICIT NONE |
---|
26 | TYPE(t_field),POINTER :: f_ps(:) |
---|
27 | TYPE(t_field),POINTER :: f_phis(:) |
---|
28 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
---|
29 | TYPE(t_field),POINTER :: f_u(:) |
---|
30 | TYPE(t_field),POINTER :: f_q(:) |
---|
31 | |
---|
32 | REAL(rstd),POINTER :: ps(:) |
---|
33 | REAL(rstd),POINTER :: phis(:) |
---|
34 | REAL(rstd),POINTER :: theta_rhodz(:,:) |
---|
35 | REAL(rstd),POINTER :: u(:,:) |
---|
36 | REAL(rstd),POINTER :: q(:,:,:) |
---|
37 | INTEGER :: ind |
---|
38 | |
---|
39 | DO ind=1,ndomain |
---|
40 | CALL swap_dimensions(ind) |
---|
41 | CALL swap_geometry(ind) |
---|
42 | ps=f_ps(ind) |
---|
43 | phis=f_phis(ind) |
---|
44 | theta_rhodz=f_theta_rhodz(ind) |
---|
45 | u=f_u(ind) |
---|
46 | q=f_q(ind) |
---|
47 | CALL compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q) |
---|
48 | ENDDO |
---|
49 | |
---|
50 | END SUBROUTINE etat0 |
---|
51 | |
---|
52 | SUBROUTINE compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q) |
---|
53 | USE icosa |
---|
54 | USE disvert_mod |
---|
55 | USE pression_mod |
---|
56 | USE exner_mod |
---|
57 | USE geopotential_mod |
---|
58 | USE theta2theta_rhodz_mod |
---|
59 | IMPLICIT NONE |
---|
60 | REAL(rstd),INTENT(OUT) :: ps(iim*jjm) |
---|
61 | REAL(rstd),INTENT(OUT) :: phis(iim*jjm) |
---|
62 | REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) |
---|
63 | REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) |
---|
64 | REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) |
---|
65 | |
---|
66 | INTEGER :: i,j,l,ij |
---|
67 | REAL(rstd) :: theta(iim*jjm,llm) |
---|
68 | REAL(rstd) :: Y(iim*jjm,llm) |
---|
69 | REAL(rstd) :: vort |
---|
70 | REAL(rstd) :: eta(llm) |
---|
71 | REAL(rstd) :: etav(llm) |
---|
72 | REAL(rstd) :: etas, etavs |
---|
73 | REAL(rstd) :: lon,lat |
---|
74 | REAL(rstd) :: ulon(3) |
---|
75 | REAL(rstd) :: ep(3), norm_ep |
---|
76 | REAL(rstd) :: Tave, T |
---|
77 | REAL(rstd) :: phis_ave |
---|
78 | REAL(rstd) :: V0(3) |
---|
79 | REAL(rstd) :: r2 |
---|
80 | REAL(rstd) :: utot |
---|
81 | REAL(rstd) :: lonx,latx |
---|
82 | REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r |
---|
83 | INTEGER :: testcase |
---|
84 | |
---|
85 | lonc=Pi/9 |
---|
86 | latc=2*Pi/9 |
---|
87 | |
---|
88 | testcase=1 |
---|
89 | CALL getin("dcmip4_testcase",testcase) |
---|
90 | |
---|
91 | |
---|
92 | DO l=1,llm |
---|
93 | eta(l)= 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) ) |
---|
94 | etav(l)=(eta(l)-eta0)*Pi/2 |
---|
95 | ENDDO |
---|
96 | etas=ap(1)/preff+bp(1) |
---|
97 | etavs=(etas-eta0)*Pi/2 |
---|
98 | |
---|
99 | DO j=jj_begin,jj_end |
---|
100 | DO i=ii_begin,ii_end |
---|
101 | ij=(j-1)*iim+i |
---|
102 | ps(ij)=ps0 |
---|
103 | ENDDO |
---|
104 | ENDDO |
---|
105 | |
---|
106 | |
---|
107 | CALL lonlat2xyz(lonc,latc,V0) |
---|
108 | |
---|
109 | u(:,:)=1e10 |
---|
110 | DO l=1,llm |
---|
111 | DO j=jj_begin-1,jj_end+1 |
---|
112 | DO i=ii_begin-1,ii_end+1 |
---|
113 | ij=(j-1)*iim+i |
---|
114 | |
---|
115 | CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) |
---|
116 | K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) |
---|
117 | r=radius*acos(K) |
---|
118 | utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) |
---|
119 | u(ij+u_right,l) = utot * sum(elon_e(ij+u_right,:) * ep_e(ij+u_right,:)) |
---|
120 | |
---|
121 | |
---|
122 | CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) |
---|
123 | K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) |
---|
124 | r=radius*acos(K) |
---|
125 | utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) |
---|
126 | u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:)) |
---|
127 | |
---|
128 | CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) |
---|
129 | K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) |
---|
130 | r=radius*acos(K) |
---|
131 | utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) |
---|
132 | u(ij+u_ldown,l) = utot * sum(elon_e(ij+u_ldown,:) * ep_e(ij+u_ldown,:)) |
---|
133 | |
---|
134 | ENDDO |
---|
135 | ENDDO |
---|
136 | ENDDO |
---|
137 | |
---|
138 | |
---|
139 | DO l=1,llm |
---|
140 | Tave=T0*eta(l)**(Rd*Gamma/g) |
---|
141 | IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 |
---|
142 | DO j=jj_begin-1,jj_end+1 |
---|
143 | DO i=ii_begin-1,ii_end+1 |
---|
144 | ij=(j-1)*iim+i |
---|
145 | CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) |
---|
146 | |
---|
147 | Y(ij,l)=((-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 & |
---|
148 | + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega) |
---|
149 | T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) |
---|
150 | |
---|
151 | theta(ij,l)=T*eta(l)**(-kappa) |
---|
152 | |
---|
153 | ENDDO |
---|
154 | ENDDO |
---|
155 | ENDDO |
---|
156 | |
---|
157 | |
---|
158 | phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) |
---|
159 | DO j=jj_begin,jj_end |
---|
160 | DO i=ii_begin,ii_end |
---|
161 | ij=(j-1)*iim+i |
---|
162 | CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) |
---|
163 | 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 & |
---|
164 | +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega ) |
---|
165 | ! phis(ij)=phis_ave+u0*cos(etavs)**1.5 |
---|
166 | |
---|
167 | ENDDO |
---|
168 | ENDDO |
---|
169 | |
---|
170 | CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) |
---|
171 | |
---|
172 | |
---|
173 | IF (testcase==1) THEN |
---|
174 | |
---|
175 | q(:,:,1)=theta(:,:) |
---|
176 | IF(nqtot>2) q(:,:,3)=1. |
---|
177 | |
---|
178 | DO l=1,llm |
---|
179 | dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1) |
---|
180 | IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa) & |
---|
181 | + kappa * (etat-eta(l))**5 * eta(l)**(-kappa-1)) |
---|
182 | DO j=jj_begin,jj_end |
---|
183 | DO i=ii_begin,ii_end |
---|
184 | ij=(j-1)*iim+i |
---|
185 | CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) |
---|
186 | dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) & |
---|
187 | + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l) & |
---|
188 | - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)& |
---|
189 | - 9./8. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l)) & |
---|
190 | * (-2*sin(lat)**6*(cos(lat)**2+1./3.)+10./63.) |
---|
191 | dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5 & |
---|
192 | *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) & |
---|
193 | + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat))) |
---|
194 | |
---|
195 | duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l)) |
---|
196 | |
---|
197 | K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) |
---|
198 | r=radius*acos(K) |
---|
199 | vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2) & |
---|
200 | + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) & |
---|
201 | -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2))) |
---|
202 | q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta)) |
---|
203 | IF(nqtot>3) q(ij,l,4)=cos(lon)*cos(lat) |
---|
204 | ENDDO |
---|
205 | ENDDO |
---|
206 | ENDDO |
---|
207 | |
---|
208 | ELSE IF (testcase==2) THEN |
---|
209 | |
---|
210 | DO l=1,llm |
---|
211 | DO j=jj_begin,jj_end |
---|
212 | DO i=ii_begin,ii_end |
---|
213 | ij=(j-1)*iim+i |
---|
214 | CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) |
---|
215 | q(ij,l,1)=q0*exp(-(lat/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2) |
---|
216 | ENDDO |
---|
217 | ENDDO |
---|
218 | ENDDO |
---|
219 | |
---|
220 | ENDIF |
---|
221 | |
---|
222 | |
---|
223 | END SUBROUTINE compute_etat0_dcmip4 |
---|
224 | |
---|
225 | END MODULE etat0_dcmip4_mod |
---|