source: codes/icosagcm/devel/src/diagnostics/wind.f90 @ 584

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

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