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

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

Simplify the management of the module.

YM
File size: 9.0 KB
Line 
1MODULE etat0_academic_mod
2
3
4
5
6CONTAINS
7 
8  SUBROUTINE test_etat0_academic
9  USE icosa
10  USE kinetic_mod
11  IMPLICIT NONE
12    TYPE(t_field),POINTER :: f_ps(:)
13    TYPE(t_field),POINTER :: f_phis(:)
14    TYPE(t_field),POINTER :: f_theta_rhodz(:)
15    TYPE(t_field),POINTER :: f_u(:)
16    TYPE(t_field),POINTER :: f_q(:)
17    TYPE(t_field),POINTER :: f_Ki(:)
18    TYPE(t_field),POINTER :: 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      CALL swap_dimensions(ind)
61      CALL swap_geometry(ind)
62      ps=f_ps(ind)
63      phis=f_phis(ind)
64      theta_rhodz=f_theta_rhodz(ind)
65      u=f_u(ind)
66      CALL compute_etat0_academic(ps, phis, theta_rhodz, u)
67    ENDDO
68
69  END SUBROUTINE etat0
70 
71  SUBROUTINE compute_etat0_academic(ps, phis, theta_rhodz, u)
72  USE icosa
73  USE disvert_mod
74  USE pression_mod
75  USE exner_mod
76  USE geopotential_mod
77  USE theta2theta_rhodz_mod
78  IMPLICIT NONE 
79  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
80  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
81  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
82  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
83 
84  INTEGER :: i,j,l,ij
85  REAL(rstd) :: r
86  REAL(rstd) :: theta(iim*jjm,llm)
87  REAL(rstd) :: zsig
88  INTEGER    :: lsup
89  REAL(rstd) :: ddsin
90  REAL(rstd) :: thetarappel
91  REAL(rstd) :: lon,lat
92  REAL(rstd) :: p(iim*jjm,llm+1)
93  REAL(rstd) :: alpha(iim*jjm,llm),beta(iim*jjm,llm)
94  REAL(rstd) :: delta
95  REAL(rstd) :: pks(iim*jjm),pk(iim*jjm,llm)
96  REAL(rstd) :: phi(iim*jjm,llm)
97  REAL(rstd) :: x 
98  REAL(rstd) :: fact(3*iim*jjm)
99  REAL(rstd) :: ut(3*iim*jjm,llm)
100 
101 
102    DO l=1,llm
103       zsig=ap(l)/preff+bp(l)
104       IF (zsig.gt.0.3) THEN
105         lsup=l
106         thetarappel=1./8.*(-log(zsig)-.5)
107         DO j=jj_begin-1,jj_end+1
108           DO i=ii_begin-1,ii_end+1
109             ij=(j-1)*iim+i
110             CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
111             ddsin=sin(lat)-sin(pi/20.)
112             theta(ij,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+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    CALL compute_pression(ps,p,1)     
137    CALL compute_exner(ps,p,pks,pk,1) 
138    CALL compute_geopotential(phis,pks,pk,theta,phi,1)
139
140
141     DO j=jj_begin-1,jj_end+1
142       DO i=ii_begin-1,ii_end+1
143         ij=(j-1)*iim+i
144
145         CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
146         IF (abs(sin(lat))<1.e-4) lat=1e-4
147         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
148         fact(ij+u_right)=(1.-x)/(2.*omega*sin(lat))
149
150         CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
151         IF (abs(sin(lat))<1.e-4) lat=1e-4
152         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
153         fact(ij+u_lup)=(1.-x)/(2.*omega*sin(lat))
154
155         CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
156         IF (abs(sin(lat))<1.e-4) lat=1e-4
157         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
158         fact(ij+u_ldown)=(1.-x)/(2.*omega*sin(lat))
159
160!         fact(ij+u_right)=-1
161!         fact(ij+u_lup)=-1
162!         fact(ij+u_ldown)=-1
163
164!         CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
165!         ps(ij)=lat
166       ENDDO
167     ENDDO
168
169 ! gradient ==> tangential component     
170     DO l=1,llm   
171       DO j=jj_begin-1,jj_end+1
172         DO i=ii_begin-1,ii_end+1
173           ij=(j-1)*iim+i
174           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) )
175           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) )
176           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) )
177         ENDDO
178       ENDDO
179     ENDDO
180     
181! attention au signe ??
182! reconstruction of perpendicular component             
183     DO l=1,llm   
184       DO j=jj_begin,jj_end
185         DO i=ii_begin,ii_end
186           ij=(j-1)*iim+i
187           u(ij+u_right,l) =  -fact(ij+u_right)/de(ij+u_right) *                                        & 
188                              ( wee(ij+u_right,1,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                        &
189                                wee(ij+u_right,2,1)*le(ij+u_lup)*ut(ij+u_lup,l)+                        &
190                                wee(ij+u_right,3,1)*le(ij+u_left)*ut(ij+u_left,l)+                      &
191                                wee(ij+u_right,4,1)*le(ij+u_ldown)*ut(ij+u_ldown,l)+                    &
192                                wee(ij+u_right,5,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                    & 
193                                wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ut(ij+t_right+u_ldown,l)+    &
194                                wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ut(ij+t_right+u_rdown,l)+    &
195                                wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ut(ij+t_right+u_right,l)+    &
196                                wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ut(ij+t_right+u_rup,l)+        &
197                                wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ut(ij+t_right+u_lup,l) )
198                               
199           u(ij+u_lup,l) = -fact(ij+u_lup)/de(ij+u_lup) *                                        & 
200                       ( wee(ij+u_lup,1,1)*le(ij+u_left)*ut(ij+u_left,l)+                        &
201                         wee(ij+u_lup,2,1)*le(ij+u_ldown)*ut(ij+u_ldown,l)+                      &
202                         wee(ij+u_lup,3,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                      &
203                         wee(ij+u_lup,4,1)*le(ij+u_right)*ut(ij+u_right,l)+                      &
204                         wee(ij+u_lup,5,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                          & 
205                         wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ut(ij+t_lup+u_right,l)+          &
206                         wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ut(ij+t_lup+u_rup,l)+              &
207                         wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ut(ij+t_lup+u_lup,l)+              &
208                         wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ut(ij+t_lup+u_left,l)+            &
209                         wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ut(ij+t_lup+u_ldown,l) )
210
211
212           u(ij+u_ldown,l) = -fact(ij+u_ldown)/de(ij+u_ldown) *                                  & 
213                       ( wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                    &
214                         wee(ij+u_ldown,2,1)*le(ij+u_right)*ut(ij+u_right,l)+                    &
215                         wee(ij+u_ldown,3,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                        &
216                         wee(ij+u_ldown,4,1)*le(ij+u_lup)*ut(ij+u_lup,l)+                        &
217                         wee(ij+u_ldown,5,1)*le(ij+u_left)*ut(ij+u_left,l)+                      & 
218                         wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ut(ij+t_ldown+u_lup,l)+        &
219                         wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ut(ij+t_ldown+u_left,l)+      &
220                         wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ut(ij+t_ldown+u_ldown,l)+    &
221                         wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ut(ij+t_ldown+u_rdown,l)+    &
222                         wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ut(ij+t_ldown+u_right,l) )
223                               
224!           u(ij+u_right,l)=fact(ij+u_right)*ut(ij+u_right,l)
225!           u(ij+u_lup,l)=fact(ij+u_lup)*ut(ij+u_lup,l)
226!           u(ij+u_ldown,l)=fact(ij+u_ldown)*ut(ij+u_ldown,l)
227!            ps(ij)=phi(ij,1)
228         ENDDO
229       ENDDO
230     ENDDO
231   
232     DO l=1,llm
233      DO j=jj_begin,jj_end
234        DO i=ii_begin,ii_end
235          ij=(j-1)*iim+i
236          CALL RANDOM_NUMBER(r) ;r=0
237          theta(ij,l)=theta(ij,l)*(1-0.0005*r)
238       ENDDO
239     ENDDO
240    ENDDO
241
242    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
243       
244
245  END SUBROUTINE compute_etat0_academic
246 
247END MODULE etat0_academic_mod
Note: See TracBrowser for help on using the repository browser.