source: codes/icosagcm/devel/src/diagnostics/wind.F90 @ 585

Last change on this file since 585 was 585, checked in by dubos, 7 years ago

devel : reconstruct fluxes at cell centers
File size: 12.4 KB
Line 
1MODULE wind_mod
2  USE omp_para
3  USE icosa
4  IMPLICIT NONE
5  PRIVATE
6
7  PUBLIC :: compute_wind_centered, compute_flux_centered, &
8       compute_wind_centered_lonlat_compound, compute_wind2d_perp_from_lonlat_compound, &
9       compute_wind_centered_from_lonlat_compound, compute_wind_perp_from_lonlat_compound, &
10       un2ulonlat, ulonlat2un
11       
12CONTAINS
13
14  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat)
15    TYPE(t_field), POINTER :: f_u(:) ! IN  : normal velocity components on edges
16    TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons   
17    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
18    INTEGER :: ind
19
20    DO ind=1,ndomain
21       IF (.NOT. assigned_domain(ind)) CYCLE
22       CALL swap_dimensions(ind)
23       CALL swap_geometry(ind)
24       u=f_u(ind)
25       ulon=f_ulon(ind)
26       ulat=f_ulat(ind)
27       CALL compute_un2ulonlat(u,ulon, ulat)
28    END DO
29
30  END SUBROUTINE un2ulonlat
31
32  SUBROUTINE ulonlat2un(f_ulon, f_ulat,f_u)
33    TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! IN : velocity reconstructed at hexagons
34    TYPE(t_field), POINTER :: f_u(:) ! OUT  : normal velocity components on edges
35   
36    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
37    INTEGER :: ind
38
39    DO ind=1,ndomain
40       IF (.NOT. assigned_domain(ind)) CYCLE
41       CALL swap_dimensions(ind)
42       CALL swap_geometry(ind)
43       u=f_u(ind)
44       ulon=f_ulon(ind)
45       ulat=f_ulat(ind)
46       CALL compute_ulonlat2un(ulon, ulat,u)
47    END DO
48
49  END SUBROUTINE ulonlat2un
50 
51  SUBROUTINE compute_wind_centered(ue,ucenter)
52  REAL(rstd) :: ue(3*iim*jjm,llm)
53  REAL(rstd) :: ucenter(iim*jjm,llm,3)
54  INTEGER :: ij,l
55  REAL(rstd) :: ue_le, cx,cy,cz, ux,uy,uz
56#include "../kernels/wind_centered.k90"
57 END SUBROUTINE compute_wind_centered
58 
59  SUBROUTINE compute_flux_centered(ue,ucenter)
60  REAL(rstd) :: ue(3*iim*jjm,llm)
61  REAL(rstd) :: ucenter(iim*jjm,llm,3)
62  INTEGER :: ij,l
63  REAL(rstd) :: ue_le, cx,cy,cz, ux,uy,uz
64#include "../kernels/flux_centered.k90"
65  END SUBROUTINE compute_flux_centered
66 
67 
68 SUBROUTINE compute_wind_on_edge(ue,uedge)
69  REAL(rstd) :: ue(3*iim*jjm,llm)
70  REAL(rstd) :: uedge(3*iim*jjm,llm,3)
71
72  REAL(rstd) :: ut(3*iim*jjm,llm)
73  INTEGER :: i,j,ij,l     
74   
75    CALL compute_tangential_compound(ue,ut)
76 
77    DO l=ll_begin,ll_end
78      DO j=jj_begin,jj_end
79        DO i=ii_begin,ii_end
80          ij=(j-1)*iim+i
81          uedge(ij+u_right,l,:)=ue(ij+u_right,l)*ep_e(ij+u_right,:)*ne(ij,right) + ut(ij+u_right,l)*et_e(ij+u_right,:)*ne(ij,right) 
82          uedge(ij+u_lup,l,:)=ue(ij+u_lup,l)*ep_e(ij+u_lup,:)*ne(ij,lup) + ut(ij+u_lup,l)*et_e(ij+u_lup,:)*ne(ij,lup)
83          uedge(ij+u_ldown,l,:)=ue(ij+u_ldown,l)*ep_e(ij+u_ldown,:)*ne(ij,ldown) + ut(ij+u_ldown,l)*et_e(ij+u_ldown,:)*ne(ij,ldown)
84        ENDDO
85      ENDDO
86    ENDDO
87 
88 END SUBROUTINE compute_wind_on_edge
89 
90 
91 
92 SUBROUTINE compute_tangential_compound(ue,ut)
93  REAL(rstd) :: ue(3*iim*jjm,llm)
94  REAL(rstd) :: ut(3*iim*jjm,llm)
95  INTEGER :: i,j,l,ij
96   
97  DO l=ll_begin,ll_end
98    DO j=jj_begin,jj_end
99      DO i=ii_begin,ii_end
100        ij=(j-1)*iim+i
101   
102        ut(ij+u_right,l) = 1/de(ij+u_right) *                                            & 
103                         ( wee(ij+u_right,1,1)*ue(ij+u_rup,l)+                           &
104                           wee(ij+u_right,2,1)*ue(ij+u_lup,l)+                           &
105                           wee(ij+u_right,3,1)*ue(ij+u_left,l)+                          &
106                           wee(ij+u_right,4,1)*ue(ij+u_ldown,l)+                         &
107                           wee(ij+u_right,5,1)*ue(ij+u_rdown,l)+                         & 
108                           wee(ij+u_right,1,2)*ue(ij+t_right+u_ldown,l)+                 &
109                           wee(ij+u_right,2,2)*ue(ij+t_right+u_rdown,l)+                 &
110                           wee(ij+u_right,3,2)*ue(ij+t_right+u_right,l)+                 &
111                           wee(ij+u_right,4,2)*ue(ij+t_right+u_rup,l)+                   &
112                           wee(ij+u_right,5,2)*ue(ij+t_right+u_lup,l) )   
113     
114        ut(ij+u_lup,l) =  1/de(ij+u_lup) *                                           & 
115                         ( wee(ij+u_lup,1,1)*ue(ij+u_left,l)+                        &
116                           wee(ij+u_lup,2,1)*ue(ij+u_ldown,l)+                       &
117                           wee(ij+u_lup,3,1)*ue(ij+u_rdown,l)+                       &
118                           wee(ij+u_lup,4,1)*ue(ij+u_right,l)+                       &
119                           wee(ij+u_lup,5,1)*ue(ij+u_rup,l)+                         & 
120                           wee(ij+u_lup,1,2)*ue(ij+t_lup+u_right,l)+                 &
121                           wee(ij+u_lup,2,2)*ue(ij+t_lup+u_rup,l)+                   &
122                           wee(ij+u_lup,3,2)*ue(ij+t_lup+u_lup,l)+                   &
123                           wee(ij+u_lup,4,2)*ue(ij+t_lup+u_left,l)+                  &
124                           wee(ij+u_lup,5,2)*ue(ij+t_lup+u_ldown,l) )
125
126   
127        ut(ij+u_ldown,l) = 1/de(ij+u_ldown) *   & 
128                         ( wee(ij+u_ldown,1,1)*ue(ij+u_rdown,l)+                      &
129                           wee(ij+u_ldown,2,1)*ue(ij+u_right,l)+                      &
130                           wee(ij+u_ldown,3,1)*ue(ij+u_rup,l)+                        &
131                           wee(ij+u_ldown,4,1)*ue(ij+u_lup,l)+                        &
132                           wee(ij+u_ldown,5,1)*ue(ij+u_left,l)+                       & 
133                           wee(ij+u_ldown,1,2)*ue(ij+t_ldown+u_lup,l)+                &
134                           wee(ij+u_ldown,2,2)*ue(ij+t_ldown+u_left,l)+               &
135                           wee(ij+u_ldown,3,2)*ue(ij+t_ldown+u_ldown,l)+              &
136                           wee(ij+u_ldown,4,2)*ue(ij+t_ldown+u_rdown,l)+              &
137                           wee(ij+u_ldown,5,2)*ue(ij+t_ldown+u_right,l) ) 
138       
139        ENDDO
140      ENDDO
141    ENDDO
142                       
143 END SUBROUTINE compute_tangential_compound
144 
145! SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat)
146!  REAL(rstd) :: u(3*iim*jjm,3,llm)
147!  REAL(rstd) :: ulon(3*iim*jjm,3,llm)
148!  REAL(rstd) :: ulat(3*iim*jjm,3,llm)
149!
150!  INTEGER :: i,j,ij,l     
151!   
152! 
153!    DO l=ll_begin,ll_end
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!          ulon(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elon_e(ij+u_right,:))*elon_e(ij+u_right,:)
158!          ulon(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elon_e(ij+u_lup,:))*elon_e(ij+u_lup,:)
159!          ulon(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elon_e(ij+u_ldown,:))*elon_e(ij+u_ldown,:)
160!         
161!          ulat(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elat_e(ij+u_right,:))*elat_e(ij+u_right,:)
162!          ulat(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elat_e(ij+u_lup,:))*elat_e(ij+u_lup,:)
163!          ulat(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elat_e(ij+u_ldown,:))*elat_e(ij+u_ldown,:)
164!         
165!        ENDDO
166!      ENDDO
167!    ENDDO
168!
169! END SUBROUTINE compute_wind_lonlat_compound
170 
171  SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u)
172  REAL(rstd) :: u(3*iim*jjm,llm,3)
173  REAL(rstd) :: ulon(3*iim*jjm,llm)
174  REAL(rstd) :: ulat(3*iim*jjm,llm)
175
176  INTEGER :: i,j,ij,l     
177 
178    DO l=ll_begin,ll_end
179      DO j=jj_begin-1,jj_end+1
180        DO i=ii_begin-1,ii_end+1
181          ij=(j-1)*iim+i
182          u(ij+u_right,l,:)=ulon(ij+u_right,l)*elon_e(ij+u_right,:)+ ulat(ij+u_right,l)*elat_e(ij+u_right,:)
183          u(ij+u_lup,l,:)=ulon(ij+u_lup,l)*elon_e(ij+u_lup,:) + ulat(ij+u_lup,l)*elat_e(ij+u_lup,:)
184          u(ij+u_ldown,l,:)=ulon(ij+u_ldown,l)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown,l)*elat_e(ij+u_ldown,:)
185        ENDDO
186      ENDDO
187    ENDDO
188 
189  END SUBROUTINE compute_wind_from_lonlat_compound
190 
191  SUBROUTINE compute_wind_centered_from_lonlat_compound(ulon, ulat, u)
192  REAL(rstd) :: u(iim*jjm,llm,3)
193  REAL(rstd) :: ulon(iim*jjm,llm)
194  REAL(rstd) :: ulat(iim*jjm,llm)
195  INTEGER :: i,j,ij,l     
196  DO l=ll_begin,ll_end
197      DO j=jj_begin-1,jj_end+1
198        DO i=ii_begin-1,ii_end+1
199          ij=(j-1)*iim+i
200          u(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ ulat(ij,l)*elat_i(ij,:)
201        ENDDO
202      ENDDO
203    ENDDO 
204  END SUBROUTINE compute_wind_centered_from_lonlat_compound
205 
206  SUBROUTINE compute_wind2D_from_lonlat_compound(ulon, ulat, u)
207  REAL(rstd) :: u(3*iim*jjm,3)
208  REAL(rstd) :: ulon(3*iim*jjm)
209  REAL(rstd) :: ulat(3*iim*jjm)
210 
211  INTEGER :: i,j,ij
212 
213  DO j=jj_begin-1,jj_end+1
214     DO i=ii_begin-1,ii_end+1
215        ij=(j-1)*iim+i
216        u(ij+u_right,:)=ulon(ij+u_right)*elon_e(ij+u_right,:)+ ulat(ij+u_right)*elat_e(ij+u_right,:)
217        u(ij+u_lup,:)=ulon(ij+u_lup)*elon_e(ij+u_lup,:) + ulat(ij+u_lup)*elat_e(ij+u_lup,:)
218        u(ij+u_ldown,:)=ulon(ij+u_ldown)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown)*elat_e(ij+u_ldown,:)
219     ENDDO
220  ENDDO
221   
222  END SUBROUTINE compute_wind2D_from_lonlat_compound
223 
224  SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up)
225  REAL(rstd) :: up(3*iim*jjm,llm)
226  REAL(rstd) :: ulon(3*iim*jjm,llm)
227  REAL(rstd) :: ulat(3*iim*jjm,llm)
228  REAL(rstd) :: u(3*iim*jjm,llm,3)
229
230  INTEGER :: i,j,ij,l     
231 
232   CALL compute_wind_from_lonlat_compound(ulon, ulat, u)
233
234    DO l=ll_begin,ll_end
235      DO j=jj_begin-1,jj_end+1
236        DO i=ii_begin-1,ii_end+1
237          ij=(j-1)*iim+i
238          up(ij+u_right,l)=sum(u(ij+u_right,l,:)*ep_e(ij+u_right,:))
239          up(ij+u_lup,l)=sum(u(ij+u_lup,l,:)*ep_e(ij+u_lup,:))
240          up(ij+u_ldown,l)=sum(u(ij+u_ldown,l,:)*ep_e(ij+u_ldown,:))
241        ENDDO
242      ENDDO
243    ENDDO
244 
245  END SUBROUTINE compute_wind_perp_from_lonlat_compound
246   
247  SUBROUTINE compute_wind2D_perp_from_lonlat_compound(ulon, ulat, up)
248  REAL(rstd) :: up(3*iim*jjm)
249  REAL(rstd) :: ulon(3*iim*jjm)
250  REAL(rstd) :: ulat(3*iim*jjm)
251  REAL(rstd) :: u(3*iim*jjm,3)
252
253  INTEGER :: i,j,ij 
254 
255  CALL compute_wind2D_from_lonlat_compound(ulon, ulat, u)
256  DO j=jj_begin-1,jj_end+1
257     DO i=ii_begin-1,ii_end+1
258        ij=(j-1)*iim+i
259        up(ij+u_right)=sum(u(ij+u_right,:)*ep_e(ij+u_right,:))
260        up(ij+u_lup)=sum(u(ij+u_lup,:)*ep_e(ij+u_lup,:))
261        up(ij+u_ldown)=sum(u(ij+u_ldown,:)*ep_e(ij+u_ldown,:))
262     ENDDO
263  ENDDO
264   
265  END SUBROUTINE compute_wind2D_perp_from_lonlat_compound
266   
267 SUBROUTINE compute_wind_centered_lonlat_compound(uc, ulon, ulat)
268  REAL(rstd) :: uc(iim*jjm,llm,3)
269  REAL(rstd) :: ulon(iim*jjm,llm)
270  REAL(rstd) :: ulat(iim*jjm,llm)
271
272  INTEGER :: i,j,ij,l     
273 
274    DO l=ll_begin,ll_end
275      DO j=jj_begin,jj_end
276        DO i=ii_begin,ii_end
277          ij=(j-1)*iim+i
278          ulon(ij,l)=sum(uc(ij,l,:)*elon_i(ij,:))
279          ulat(ij,l)=sum(uc(ij,l,:)*elat_i(ij,:)) 
280        ENDDO
281      ENDDO
282    ENDDO
283 
284 END SUBROUTINE compute_wind_centered_lonlat_compound
285
286 SUBROUTINE compute_wind_centered_from_wind_lonlat_centered(ulon, ulat,uc)
287  REAL(rstd) :: ulon(iim*jjm,llm)
288  REAL(rstd) :: ulat(iim*jjm,llm)
289  REAL(rstd) :: uc(iim*jjm,llm,3)
290
291  INTEGER :: i,j,ij,l     
292   
293 
294    DO l=ll_begin,ll_end
295      DO j=jj_begin,jj_end
296        DO i=ii_begin,ii_end
297          ij=(j-1)*iim+i
298          uc(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ulat(ij,l)*elat_i(ij,:)
299        ENDDO
300      ENDDO
301    ENDDO
302 
303 END SUBROUTINE compute_wind_centered_from_wind_lonlat_centered
304
305 SUBROUTINE compute_wind_perp_from_wind_centered(uc,un)
306
307  IMPLICIT NONE
308  REAL(rstd),INTENT(IN)   :: uc(iim*jjm,llm,3)
309  REAL(rstd),INTENT(OUT)  :: un(3*iim*jjm,llm)
310
311  INTEGER :: i,j,ij,l     
312   
313 
314    DO l=ll_begin,ll_end
315      DO j=jj_begin,jj_end
316        DO i=ii_begin,ii_end
317          ij=(j-1)*iim+i
318          un(ij+u_right,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_right,l,:))*ep_e(ij+u_right,:))
319          un(ij+u_lup,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_lup,l,:))*ep_e(ij+u_lup,:))
320          un(ij+u_ldown,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:))
321         ENDDO
322      ENDDO
323    ENDDO
324 
325 END SUBROUTINE compute_wind_perp_from_wind_centered
326
327
328 SUBROUTINE compute_un2ulonlat(un, ulon, ulat)
329  REAL(rstd),INTENT(IN)  :: un(3*iim*jjm,llm)
330  REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm)
331  REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm)
332
333  REAL(rstd)             :: uc(iim*jjm,llm,3)
334   
335  CALL compute_wind_centered(un,uc) 
336  CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
337 
338 END SUBROUTINE compute_un2ulonlat
339
340 SUBROUTINE compute_ulonlat2un(ulon, ulat,un)
341  REAL(rstd),INTENT(IN) :: ulon(iim*jjm,llm)
342  REAL(rstd),INTENT(IN) :: ulat(iim*jjm,llm)
343  REAL(rstd),INTENT(OUT)  :: un(3*iim*jjm,llm)
344
345  REAL(rstd)             :: uc(iim*jjm,llm,3)
346   
347    CALL compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc) 
348    CALL compute_wind_perp_from_wind_centered(uc, un)
349 
350 END SUBROUTINE compute_ulonlat2un
351
352
353END MODULE wind_mod
Note: See TracBrowser for help on using the repository browser.