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

Last change on this file since 238 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

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,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(3*iim*jjm,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(iim*jjm,llm),beta(iim*jjm,llm)
95  REAL(rstd) :: delta
96  REAL(rstd) :: pks(iim*jjm),pk(iim*jjm,llm)
97  REAL(rstd) :: phi(iim*jjm,llm)
98  REAL(rstd) :: x 
99  REAL(rstd) :: fact(3*iim*jjm)
100  REAL(rstd) :: ut(3*iim*jjm,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             CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
112             ddsin=sin(lat)-sin(pi/20.)
113             theta(ij,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+thetarappel)
114            ENDDO
115         ENDDO
116       ELSE
117!   Choix isotherme au-dessus de 300 mbar
118         DO j=jj_begin-1,jj_end+1
119           DO i=ii_begin-1,ii_end+1
120             ij=(j-1)*iim+i
121             theta(ij,l)=theta(ij,lsup)*(0.3/zsig)**kappa
122           ENDDO
123         ENDDO
124       ENDIF
125     ENDDO
126     
127     
128     DO j=jj_begin-1,jj_end+1
129       DO i=ii_begin-1,ii_end+1
130         ij=(j-1)*iim+i
131         ps(ij)=1.e5
132         phis(ij)=0
133       ENDDO
134     ENDDO
135     
136
137    CALL compute_pression(ps,p,1)     
138    CALL compute_exner(ps,p,pks,pk,1) 
139    CALL compute_geopotential(phis,pks,pk,theta,phi,1)
140
141
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
146         CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
147         IF (abs(sin(lat))<1.e-4) lat=1e-4
148         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
149         fact(ij+u_right)=(1.-x)/(2.*omega*sin(lat))
150
151         CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
152         IF (abs(sin(lat))<1.e-4) lat=1e-4
153         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
154         fact(ij+u_lup)=(1.-x)/(2.*omega*sin(lat))
155
156         CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
157         IF (abs(sin(lat))<1.e-4) lat=1e-4
158         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
159         fact(ij+u_ldown)=(1.-x)/(2.*omega*sin(lat))
160
161!         fact(ij+u_right)=-1
162!         fact(ij+u_lup)=-1
163!         fact(ij+u_ldown)=-1
164
165!         CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
166!         ps(ij)=lat
167       ENDDO
168     ENDDO
169
170 ! gradient ==> tangential component     
171     DO l=1,llm   
172       DO j=jj_begin-1,jj_end+1
173         DO i=ii_begin-1,ii_end+1
174           ij=(j-1)*iim+i
175           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) )
176           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) )
177           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) )
178         ENDDO
179       ENDDO
180     ENDDO
181     
182! attention au signe ??
183! reconstruction of perpendicular component             
184     DO l=1,llm   
185       DO j=jj_begin,jj_end
186         DO i=ii_begin,ii_end
187           ij=(j-1)*iim+i
188           u(ij+u_right,l) =  -fact(ij+u_right)/de(ij+u_right) *                                        & 
189                              ( wee(ij+u_right,1,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                        &
190                                wee(ij+u_right,2,1)*le(ij+u_lup)*ut(ij+u_lup,l)+                        &
191                                wee(ij+u_right,3,1)*le(ij+u_left)*ut(ij+u_left,l)+                      &
192                                wee(ij+u_right,4,1)*le(ij+u_ldown)*ut(ij+u_ldown,l)+                    &
193                                wee(ij+u_right,5,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                    & 
194                                wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ut(ij+t_right+u_ldown,l)+    &
195                                wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ut(ij+t_right+u_rdown,l)+    &
196                                wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ut(ij+t_right+u_right,l)+    &
197                                wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ut(ij+t_right+u_rup,l)+        &
198                                wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ut(ij+t_right+u_lup,l) )
199                               
200           u(ij+u_lup,l) = -fact(ij+u_lup)/de(ij+u_lup) *                                        & 
201                       ( wee(ij+u_lup,1,1)*le(ij+u_left)*ut(ij+u_left,l)+                        &
202                         wee(ij+u_lup,2,1)*le(ij+u_ldown)*ut(ij+u_ldown,l)+                      &
203                         wee(ij+u_lup,3,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                      &
204                         wee(ij+u_lup,4,1)*le(ij+u_right)*ut(ij+u_right,l)+                      &
205                         wee(ij+u_lup,5,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                          & 
206                         wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ut(ij+t_lup+u_right,l)+          &
207                         wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ut(ij+t_lup+u_rup,l)+              &
208                         wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ut(ij+t_lup+u_lup,l)+              &
209                         wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ut(ij+t_lup+u_left,l)+            &
210                         wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ut(ij+t_lup+u_ldown,l) )
211
212
213           u(ij+u_ldown,l) = -fact(ij+u_ldown)/de(ij+u_ldown) *                                  & 
214                       ( wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                    &
215                         wee(ij+u_ldown,2,1)*le(ij+u_right)*ut(ij+u_right,l)+                    &
216                         wee(ij+u_ldown,3,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                        &
217                         wee(ij+u_ldown,4,1)*le(ij+u_lup)*ut(ij+u_lup,l)+                        &
218                         wee(ij+u_ldown,5,1)*le(ij+u_left)*ut(ij+u_left,l)+                      & 
219                         wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ut(ij+t_ldown+u_lup,l)+        &
220                         wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ut(ij+t_ldown+u_left,l)+      &
221                         wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ut(ij+t_ldown+u_ldown,l)+    &
222                         wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ut(ij+t_ldown+u_rdown,l)+    &
223                         wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ut(ij+t_ldown+u_right,l) )
224                               
225!           u(ij+u_right,l)=fact(ij+u_right)*ut(ij+u_right,l)
226!           u(ij+u_lup,l)=fact(ij+u_lup)*ut(ij+u_lup,l)
227!           u(ij+u_ldown,l)=fact(ij+u_ldown)*ut(ij+u_ldown,l)
228!            ps(ij)=phi(ij,1)
229         ENDDO
230       ENDDO
231     ENDDO
232   
233     DO l=1,llm
234      DO j=jj_begin,jj_end
235        DO i=ii_begin,ii_end
236          ij=(j-1)*iim+i
237          CALL RANDOM_NUMBER(r) ;r=0
238          theta(ij,l)=theta(ij,l)*(1-0.0005*r)
239       ENDDO
240     ENDDO
241    ENDDO
242
243    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
244       
245
246  END SUBROUTINE compute_etat0_academic
247 
248END MODULE etat0_academic_mod
Note: See TracBrowser for help on using the repository browser.