Context Navigation

etat0_academic.f90 @ 3850

Last change on this file since 3850 was 3810, checked in by ymipsl, 11 years ago
Add DYNAMICO in aquaplanet configuration YM
File size: 8.8 KB

Line
1	MODULE etat0_academic_mod
2
3
4
5
6	CONTAINS
7
8	SUBROUTINE test_etat0_academic
9	USE icosa
10	USE kinetic_mod
11	IMPLICIT NONE
12	TYPE(t_field),POINTER,SAVE :: f_ps(:)
13	TYPE(t_field),POINTER,SAVE :: f_phis(:)
14	TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:)
15	TYPE(t_field),POINTER,SAVE :: f_u(:)
16	TYPE(t_field),POINTER,SAVE :: f_q(:)
17	TYPE(t_field),POINTER,SAVE :: f_Ki(:)
18	TYPE(t_field),POINTER,SAVE :: f_temp(:)
19
20	REAL(rstd),POINTER :: Ki(:,:)
21	REAL(rstd),POINTER :: temp(:)
22	INTEGER :: ind
23
24
25	CALL allocate_field(f_ps,field_t,type_real)
26	CALL allocate_field(f_phis,field_t,type_real)
27	CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
28	CALL allocate_field(f_u,field_u,type_real,llm)
29	CALL allocate_field(f_Ki,field_t,type_real,llm)
30	CALL allocate_field(f_temp,field_t,type_real)
31
32	CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
33
34	CALL kinetic(f_u,f_Ki)
35
36	CALL writefield('ps',f_ps)
37	CALL writefield('theta',f_theta_rhodz)
38
39	CALL writefield('Ki',f_Ki)
40
41	END SUBROUTINE test_etat0_academic
42
43
44	SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
45	USE icosa
46	IMPLICIT NONE
47	TYPE(t_field),POINTER :: f_ps(:)
48	TYPE(t_field),POINTER :: f_phis(:)
49	TYPE(t_field),POINTER :: f_theta_rhodz(:)
50	TYPE(t_field),POINTER :: f_u(:)
51	TYPE(t_field),POINTER :: f_q(:)
52
53	REAL(rstd),POINTER :: ps(:)
54	REAL(rstd),POINTER :: phis(:)
55	REAL(rstd),POINTER :: theta_rhodz(:,:)
56	REAL(rstd),POINTER :: u(:,:)
57	INTEGER :: ind
58
59	DO ind=1,ndomain
60	IF (.NOT. assigned_domain(ind)) CYCLE
61	CALL swap_dimensions(ind)
62	CALL swap_geometry(ind)
63	ps=f_ps(ind)
64	phis=f_phis(ind)
65	theta_rhodz=f_theta_rhodz(ind)
66	u=f_u(ind)
67	CALL compute_etat0_academic(ps, phis, theta_rhodz, u)
68	ENDDO
69
70	END SUBROUTINE etat0
71
72	SUBROUTINE compute_etat0_academic(ps, phis, theta_rhodz, u)
73	USE icosa
74	USE disvert_mod
75	USE pression_mod
76	USE exner_mod
77	USE geopotential_mod
78	USE theta2theta_rhodz_mod
79	IMPLICIT NONE
80	REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
81	REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
82	REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
83	REAL(rstd),INTENT(OUT) :: u(3iimjjm,llm)
84
85	INTEGER :: i,j,l,ij
86	REAL(rstd) :: r
87	REAL(rstd) :: theta(iim*jjm,llm)
88	REAL(rstd) :: zsig
89	INTEGER :: lsup
90	REAL(rstd) :: ddsin
91	REAL(rstd) :: thetarappel
92	REAL(rstd) :: lon,lat
93	REAL(rstd) :: p(iim*jjm,llm+1)
94	REAL(rstd) :: alpha(iimjjm,llm),beta(iimjjm,llm)
95	REAL(rstd) :: delta
96	REAL(rstd) :: pks(iimjjm),pk(iimjjm,llm)
97	REAL(rstd) :: phi(iim*jjm,llm)
98	REAL(rstd) :: x
99	REAL(rstd) :: fact(3iimjjm)
100	REAL(rstd) :: ut(3iimjjm,llm)
101
102
103	DO l=1,llm
104	zsig=ap(l)/preff+bp(l)
105	IF (zsig.gt.0.3) THEN
106	lsup=l
107	thetarappel=1./8.*(-log(zsig)-.5)
108	DO j=jj_begin-1,jj_end+1
109	DO i=ii_begin-1,ii_end+1
110	ij=(j-1)*iim+i
111	ddsin=sin(lat_i(ij))-sin(pi/20.)
112	theta(ij,l)=300.(1+1./18.(1.-3.ddsinddsin)+thetarappel)
113	ENDDO
114	ENDDO
115	ELSE
116	! Choix isotherme au-dessus de 300 mbar
117	DO j=jj_begin-1,jj_end+1
118	DO i=ii_begin-1,ii_end+1
119	ij=(j-1)*iim+i
120	theta(ij,l)=theta(ij,lsup)(0.3/zsig)*kappa
121	ENDDO
122	ENDDO
123	ENDIF
124	ENDDO
125
126
127	DO j=jj_begin-1,jj_end+1
128	DO i=ii_begin-1,ii_end+1
129	ij=(j-1)*iim+i
130	ps(ij)=1.e5
131	phis(ij)=0
132	ENDDO
133	ENDDO
134
135
136	!$OMP BARRIER
137	CALL compute_pression(ps,p,1)
138	!$OMP BARRIER
139	CALL compute_exner(ps,p,pks,pk,1)
140	!$OMP BARRIER
141	CALL compute_geopotential(phis,pks,pk,theta,phi,1)
142
143
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
148	lat=lat_e(ij+u_right)
149	IF (abs(sin(lat))<1.e-4) lat=1e-4
150	x=cos(lat) ; x=xx ; x=xx ; x=x*x
151	fact(ij+u_right)=(1.-x)/(2.omegasin(lat))
152
153	lat=lat_e(ij+u_lup)
154	IF (abs(sin(lat))<1.e-4) lat=1e-4
155	x=cos(lat) ; x=xx ; x=xx ; x=x*x
156	fact(ij+u_lup)=(1.-x)/(2.omegasin(lat))
157
158	lat=lat_e(ij+u_ldown)
159	IF (abs(sin(lat))<1.e-4) lat=1e-4
160	x=cos(lat) ; x=xx ; x=xx ; x=x*x
161	fact(ij+u_ldown)=(1.-x)/(2.omegasin(lat))
162	ENDDO
163	ENDDO
164
165	! gradient ==> tangential component
166	DO l=1,llm
167	DO j=jj_begin-1,jj_end+1
168	DO i=ii_begin-1,ii_end+1
169	ij=(j-1)*iim+i
170	ut(ij+u_right,l)=1/de(ij+u_right)(ne(ij,right)phi(ij,l)+ ne(ij+t_right,left)*phi(ij+t_right,l) )
171	ut(ij+u_lup,l)=1/de(ij+u_lup)(ne(ij,lup)phi(ij,l)+ ne(ij+t_lup,rdown)*phi(ij+t_lup,l) )
172	ut(ij+u_ldown,l)=1/de(ij+u_ldown)(ne(ij,ldown)phi(ij,l)+ ne(ij+t_ldown,rup)*phi(ij+t_ldown,l) )
173	ENDDO
174	ENDDO
175	ENDDO
176
177	! attention au signe ??
178	! reconstruction of perpendicular component
179	DO l=1,llm
180	DO j=jj_begin,jj_end
181	DO i=ii_begin,ii_end
182	ij=(j-1)*iim+i
183	u(ij+u_right,l) = -fact(ij+u_right)/de(ij+u_right) * &
184	( wee(ij+u_right,1,1)le(ij+u_rup)ut(ij+u_rup,l)+ &
185	wee(ij+u_right,2,1)le(ij+u_lup)ut(ij+u_lup,l)+ &
186	wee(ij+u_right,3,1)le(ij+u_left)ut(ij+u_left,l)+ &
187	wee(ij+u_right,4,1)le(ij+u_ldown)ut(ij+u_ldown,l)+ &
188	wee(ij+u_right,5,1)le(ij+u_rdown)ut(ij+u_rdown,l)+ &
189	wee(ij+u_right,1,2)le(ij+t_right+u_ldown)ut(ij+t_right+u_ldown,l)+ &
190	wee(ij+u_right,2,2)le(ij+t_right+u_rdown)ut(ij+t_right+u_rdown,l)+ &
191	wee(ij+u_right,3,2)le(ij+t_right+u_right)ut(ij+t_right+u_right,l)+ &
192	wee(ij+u_right,4,2)le(ij+t_right+u_rup)ut(ij+t_right+u_rup,l)+ &
193	wee(ij+u_right,5,2)le(ij+t_right+u_lup)ut(ij+t_right+u_lup,l) )
194
195	u(ij+u_lup,l) = -fact(ij+u_lup)/de(ij+u_lup) * &
196	( wee(ij+u_lup,1,1)le(ij+u_left)ut(ij+u_left,l)+ &
197	wee(ij+u_lup,2,1)le(ij+u_ldown)ut(ij+u_ldown,l)+ &
198	wee(ij+u_lup,3,1)le(ij+u_rdown)ut(ij+u_rdown,l)+ &
199	wee(ij+u_lup,4,1)le(ij+u_right)ut(ij+u_right,l)+ &
200	wee(ij+u_lup,5,1)le(ij+u_rup)ut(ij+u_rup,l)+ &
201	wee(ij+u_lup,1,2)le(ij+t_lup+u_right)ut(ij+t_lup+u_right,l)+ &
202	wee(ij+u_lup,2,2)le(ij+t_lup+u_rup)ut(ij+t_lup+u_rup,l)+ &
203	wee(ij+u_lup,3,2)le(ij+t_lup+u_lup)ut(ij+t_lup+u_lup,l)+ &
204	wee(ij+u_lup,4,2)le(ij+t_lup+u_left)ut(ij+t_lup+u_left,l)+ &
205	wee(ij+u_lup,5,2)le(ij+t_lup+u_ldown)ut(ij+t_lup+u_ldown,l) )
206
207
208	u(ij+u_ldown,l) = -fact(ij+u_ldown)/de(ij+u_ldown) * &
209	( wee(ij+u_ldown,1,1)le(ij+u_rdown)ut(ij+u_rdown,l)+ &
210	wee(ij+u_ldown,2,1)le(ij+u_right)ut(ij+u_right,l)+ &
211	wee(ij+u_ldown,3,1)le(ij+u_rup)ut(ij+u_rup,l)+ &
212	wee(ij+u_ldown,4,1)le(ij+u_lup)ut(ij+u_lup,l)+ &
213	wee(ij+u_ldown,5,1)le(ij+u_left)ut(ij+u_left,l)+ &
214	wee(ij+u_ldown,1,2)le(ij+t_ldown+u_lup)ut(ij+t_ldown+u_lup,l)+ &
215	wee(ij+u_ldown,2,2)le(ij+t_ldown+u_left)ut(ij+t_ldown+u_left,l)+ &
216	wee(ij+u_ldown,3,2)le(ij+t_ldown+u_ldown)ut(ij+t_ldown+u_ldown,l)+ &
217	wee(ij+u_ldown,4,2)le(ij+t_ldown+u_rdown)ut(ij+t_ldown+u_rdown,l)+ &
218	wee(ij+u_ldown,5,2)le(ij+t_ldown+u_right)ut(ij+t_ldown+u_right,l) )
219
220	! u(ij+u_right,l)=fact(ij+u_right)*ut(ij+u_right,l)
221	! u(ij+u_lup,l)=fact(ij+u_lup)*ut(ij+u_lup,l)
222	! u(ij+u_ldown,l)=fact(ij+u_ldown)*ut(ij+u_ldown,l)
223	! ps(ij)=phi(ij,1)
224	ENDDO
225	ENDDO
226	ENDDO
227
228	DO l=1,llm
229	DO j=jj_begin,jj_end
230	DO i=ii_begin,ii_end
231	ij=(j-1)*iim+i
232	CALL RANDOM_NUMBER(r) ;r=0
233	theta(ij,l)=theta(ij,l)(1-0.0005r)
234	ENDDO
235	ENDDO
236	ENDDO
237
238	CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
239
240
241	END SUBROUTINE compute_etat0_academic
242
243	END MODULE etat0_academic_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format