source: codes/icosagcm/trunk/src/etat0_academic.f90 @ 13

Last change on this file since 13 was 13, checked in by ymipsl, 12 years ago

some update

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