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

Last change on this file since 612 was 612, checked in by dubos, 6 years ago

devel : rename directory kernels as kernels_hex
File size: 12.5 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), PARAMETER :: scale=1.
56  REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz
57#include "../kernels_hex/wind_centered.k90"
58 END SUBROUTINE compute_wind_centered
59 
60  SUBROUTINE compute_flux_centered(scale,ue,ucenter)
61  REAL(rstd), INTENT(IN) :: scale
62  REAL(rstd) :: ue(3*iim*jjm,llm)
63  REAL(rstd) :: ucenter(iim*jjm,llm,3)
64  INTEGER :: ij,l
65  REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz
66#include "../kernels_hex/flux_centered.k90"
67  END SUBROUTINE compute_flux_centered
68 
69 
70 SUBROUTINE compute_wind_on_edge(ue,uedge)
71  REAL(rstd) :: ue(3*iim*jjm,llm)
72  REAL(rstd) :: uedge(3*iim*jjm,llm,3)
73
74  REAL(rstd) :: ut(3*iim*jjm,llm)
75  INTEGER :: i,j,ij,l     
76   
77    CALL compute_tangential_compound(ue,ut)
78 
79    DO l=ll_begin,ll_end
80      DO j=jj_begin,jj_end
81        DO i=ii_begin,ii_end
82          ij=(j-1)*iim+i
83          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) 
84          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)
85          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)
86        ENDDO
87      ENDDO
88    ENDDO
89 
90 END SUBROUTINE compute_wind_on_edge
91 
92 
93 
94 SUBROUTINE compute_tangential_compound(ue,ut)
95  REAL(rstd) :: ue(3*iim*jjm,llm)
96  REAL(rstd) :: ut(3*iim*jjm,llm)
97  INTEGER :: i,j,l,ij
98   
99  DO l=ll_begin,ll_end
100    DO j=jj_begin,jj_end
101      DO i=ii_begin,ii_end
102        ij=(j-1)*iim+i
103   
104        ut(ij+u_right,l) = 1/de(ij+u_right) *                                            & 
105                         ( wee(ij+u_right,1,1)*ue(ij+u_rup,l)+                           &
106                           wee(ij+u_right,2,1)*ue(ij+u_lup,l)+                           &
107                           wee(ij+u_right,3,1)*ue(ij+u_left,l)+                          &
108                           wee(ij+u_right,4,1)*ue(ij+u_ldown,l)+                         &
109                           wee(ij+u_right,5,1)*ue(ij+u_rdown,l)+                         & 
110                           wee(ij+u_right,1,2)*ue(ij+t_right+u_ldown,l)+                 &
111                           wee(ij+u_right,2,2)*ue(ij+t_right+u_rdown,l)+                 &
112                           wee(ij+u_right,3,2)*ue(ij+t_right+u_right,l)+                 &
113                           wee(ij+u_right,4,2)*ue(ij+t_right+u_rup,l)+                   &
114                           wee(ij+u_right,5,2)*ue(ij+t_right+u_lup,l) )   
115     
116        ut(ij+u_lup,l) =  1/de(ij+u_lup) *                                           & 
117                         ( wee(ij+u_lup,1,1)*ue(ij+u_left,l)+                        &
118                           wee(ij+u_lup,2,1)*ue(ij+u_ldown,l)+                       &
119                           wee(ij+u_lup,3,1)*ue(ij+u_rdown,l)+                       &
120                           wee(ij+u_lup,4,1)*ue(ij+u_right,l)+                       &
121                           wee(ij+u_lup,5,1)*ue(ij+u_rup,l)+                         & 
122                           wee(ij+u_lup,1,2)*ue(ij+t_lup+u_right,l)+                 &
123                           wee(ij+u_lup,2,2)*ue(ij+t_lup+u_rup,l)+                   &
124                           wee(ij+u_lup,3,2)*ue(ij+t_lup+u_lup,l)+                   &
125                           wee(ij+u_lup,4,2)*ue(ij+t_lup+u_left,l)+                  &
126                           wee(ij+u_lup,5,2)*ue(ij+t_lup+u_ldown,l) )
127
128   
129        ut(ij+u_ldown,l) = 1/de(ij+u_ldown) *   & 
130                         ( wee(ij+u_ldown,1,1)*ue(ij+u_rdown,l)+                      &
131                           wee(ij+u_ldown,2,1)*ue(ij+u_right,l)+                      &
132                           wee(ij+u_ldown,3,1)*ue(ij+u_rup,l)+                        &
133                           wee(ij+u_ldown,4,1)*ue(ij+u_lup,l)+                        &
134                           wee(ij+u_ldown,5,1)*ue(ij+u_left,l)+                       & 
135                           wee(ij+u_ldown,1,2)*ue(ij+t_ldown+u_lup,l)+                &
136                           wee(ij+u_ldown,2,2)*ue(ij+t_ldown+u_left,l)+               &
137                           wee(ij+u_ldown,3,2)*ue(ij+t_ldown+u_ldown,l)+              &
138                           wee(ij+u_ldown,4,2)*ue(ij+t_ldown+u_rdown,l)+              &
139                           wee(ij+u_ldown,5,2)*ue(ij+t_ldown+u_right,l) ) 
140       
141        ENDDO
142      ENDDO
143    ENDDO
144                       
145 END SUBROUTINE compute_tangential_compound
146 
147! SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat)
148!  REAL(rstd) :: u(3*iim*jjm,3,llm)
149!  REAL(rstd) :: ulon(3*iim*jjm,3,llm)
150!  REAL(rstd) :: ulat(3*iim*jjm,3,llm)
151!
152!  INTEGER :: i,j,ij,l     
153!   
154! 
155!    DO l=ll_begin,ll_end
156!      DO j=jj_begin-1,jj_end+1
157!        DO i=ii_begin-1,ii_end+1
158!          ij=(j-1)*iim+i
159!          ulon(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elon_e(ij+u_right,:))*elon_e(ij+u_right,:)
160!          ulon(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elon_e(ij+u_lup,:))*elon_e(ij+u_lup,:)
161!          ulon(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elon_e(ij+u_ldown,:))*elon_e(ij+u_ldown,:)
162!         
163!          ulat(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elat_e(ij+u_right,:))*elat_e(ij+u_right,:)
164!          ulat(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elat_e(ij+u_lup,:))*elat_e(ij+u_lup,:)
165!          ulat(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elat_e(ij+u_ldown,:))*elat_e(ij+u_ldown,:)
166!         
167!        ENDDO
168!      ENDDO
169!    ENDDO
170!
171! END SUBROUTINE compute_wind_lonlat_compound
172 
173  SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u)
174  REAL(rstd) :: u(3*iim*jjm,llm,3)
175  REAL(rstd) :: ulon(3*iim*jjm,llm)
176  REAL(rstd) :: ulat(3*iim*jjm,llm)
177
178  INTEGER :: i,j,ij,l     
179 
180    DO l=ll_begin,ll_end
181      DO j=jj_begin-1,jj_end+1
182        DO i=ii_begin-1,ii_end+1
183          ij=(j-1)*iim+i
184          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,:)
185          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,:)
186          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,:)
187        ENDDO
188      ENDDO
189    ENDDO
190 
191  END SUBROUTINE compute_wind_from_lonlat_compound
192 
193  SUBROUTINE compute_wind_centered_from_lonlat_compound(ulon, ulat, u)
194  REAL(rstd) :: u(iim*jjm,llm,3)
195  REAL(rstd) :: ulon(iim*jjm,llm)
196  REAL(rstd) :: ulat(iim*jjm,llm)
197  INTEGER :: i,j,ij,l     
198  DO l=ll_begin,ll_end
199      DO j=jj_begin-1,jj_end+1
200        DO i=ii_begin-1,ii_end+1
201          ij=(j-1)*iim+i
202          u(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ ulat(ij,l)*elat_i(ij,:)
203        ENDDO
204      ENDDO
205    ENDDO 
206  END SUBROUTINE compute_wind_centered_from_lonlat_compound
207 
208  SUBROUTINE compute_wind2D_from_lonlat_compound(ulon, ulat, u)
209  REAL(rstd) :: u(3*iim*jjm,3)
210  REAL(rstd) :: ulon(3*iim*jjm)
211  REAL(rstd) :: ulat(3*iim*jjm)
212 
213  INTEGER :: i,j,ij
214 
215  DO j=jj_begin-1,jj_end+1
216     DO i=ii_begin-1,ii_end+1
217        ij=(j-1)*iim+i
218        u(ij+u_right,:)=ulon(ij+u_right)*elon_e(ij+u_right,:)+ ulat(ij+u_right)*elat_e(ij+u_right,:)
219        u(ij+u_lup,:)=ulon(ij+u_lup)*elon_e(ij+u_lup,:) + ulat(ij+u_lup)*elat_e(ij+u_lup,:)
220        u(ij+u_ldown,:)=ulon(ij+u_ldown)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown)*elat_e(ij+u_ldown,:)
221     ENDDO
222  ENDDO
223   
224  END SUBROUTINE compute_wind2D_from_lonlat_compound
225 
226  SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up)
227  REAL(rstd) :: up(3*iim*jjm,llm)
228  REAL(rstd) :: ulon(3*iim*jjm,llm)
229  REAL(rstd) :: ulat(3*iim*jjm,llm)
230  REAL(rstd) :: u(3*iim*jjm,llm,3)
231
232  INTEGER :: i,j,ij,l     
233 
234   CALL compute_wind_from_lonlat_compound(ulon, ulat, u)
235
236    DO l=ll_begin,ll_end
237      DO j=jj_begin-1,jj_end+1
238        DO i=ii_begin-1,ii_end+1
239          ij=(j-1)*iim+i
240          up(ij+u_right,l)=sum(u(ij+u_right,l,:)*ep_e(ij+u_right,:))
241          up(ij+u_lup,l)=sum(u(ij+u_lup,l,:)*ep_e(ij+u_lup,:))
242          up(ij+u_ldown,l)=sum(u(ij+u_ldown,l,:)*ep_e(ij+u_ldown,:))
243        ENDDO
244      ENDDO
245    ENDDO
246 
247  END SUBROUTINE compute_wind_perp_from_lonlat_compound
248   
249  SUBROUTINE compute_wind2D_perp_from_lonlat_compound(ulon, ulat, up)
250  REAL(rstd) :: up(3*iim*jjm)
251  REAL(rstd) :: ulon(3*iim*jjm)
252  REAL(rstd) :: ulat(3*iim*jjm)
253  REAL(rstd) :: u(3*iim*jjm,3)
254
255  INTEGER :: i,j,ij 
256 
257  CALL compute_wind2D_from_lonlat_compound(ulon, ulat, u)
258  DO j=jj_begin-1,jj_end+1
259     DO i=ii_begin-1,ii_end+1
260        ij=(j-1)*iim+i
261        up(ij+u_right)=sum(u(ij+u_right,:)*ep_e(ij+u_right,:))
262        up(ij+u_lup)=sum(u(ij+u_lup,:)*ep_e(ij+u_lup,:))
263        up(ij+u_ldown)=sum(u(ij+u_ldown,:)*ep_e(ij+u_ldown,:))
264     ENDDO
265  ENDDO
266   
267  END SUBROUTINE compute_wind2D_perp_from_lonlat_compound
268   
269 SUBROUTINE compute_wind_centered_lonlat_compound(uc, ulon, ulat)
270  REAL(rstd) :: uc(iim*jjm,llm,3)
271  REAL(rstd) :: ulon(iim*jjm,llm)
272  REAL(rstd) :: ulat(iim*jjm,llm)
273
274  INTEGER :: i,j,ij,l     
275 
276    DO l=ll_begin,ll_end
277      DO j=jj_begin,jj_end
278        DO i=ii_begin,ii_end
279          ij=(j-1)*iim+i
280          ulon(ij,l)=sum(uc(ij,l,:)*elon_i(ij,:))
281          ulat(ij,l)=sum(uc(ij,l,:)*elat_i(ij,:)) 
282        ENDDO
283      ENDDO
284    ENDDO
285 
286 END SUBROUTINE compute_wind_centered_lonlat_compound
287
288 SUBROUTINE compute_wind_centered_from_wind_lonlat_centered(ulon, ulat,uc)
289  REAL(rstd) :: ulon(iim*jjm,llm)
290  REAL(rstd) :: ulat(iim*jjm,llm)
291  REAL(rstd) :: uc(iim*jjm,llm,3)
292
293  INTEGER :: i,j,ij,l     
294   
295 
296    DO l=ll_begin,ll_end
297      DO j=jj_begin,jj_end
298        DO i=ii_begin,ii_end
299          ij=(j-1)*iim+i
300          uc(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ulat(ij,l)*elat_i(ij,:)
301        ENDDO
302      ENDDO
303    ENDDO
304 
305 END SUBROUTINE compute_wind_centered_from_wind_lonlat_centered
306
307 SUBROUTINE compute_wind_perp_from_wind_centered(uc,un)
308
309  IMPLICIT NONE
310  REAL(rstd),INTENT(IN)   :: uc(iim*jjm,llm,3)
311  REAL(rstd),INTENT(OUT)  :: un(3*iim*jjm,llm)
312
313  INTEGER :: i,j,ij,l     
314   
315 
316    DO l=ll_begin,ll_end
317      DO j=jj_begin,jj_end
318        DO i=ii_begin,ii_end
319          ij=(j-1)*iim+i
320          un(ij+u_right,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_right,l,:))*ep_e(ij+u_right,:))
321          un(ij+u_lup,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_lup,l,:))*ep_e(ij+u_lup,:))
322          un(ij+u_ldown,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:))
323         ENDDO
324      ENDDO
325    ENDDO
326 
327 END SUBROUTINE compute_wind_perp_from_wind_centered
328
329
330 SUBROUTINE compute_un2ulonlat(un, ulon, ulat)
331  REAL(rstd),INTENT(IN)  :: un(3*iim*jjm,llm)
332  REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm)
333  REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm)
334
335  REAL(rstd)             :: uc(iim*jjm,llm,3)
336   
337  CALL compute_wind_centered(un,uc) 
338  CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
339 
340 END SUBROUTINE compute_un2ulonlat
341
342 SUBROUTINE compute_ulonlat2un(ulon, ulat,un)
343  REAL(rstd),INTENT(IN) :: ulon(iim*jjm,llm)
344  REAL(rstd),INTENT(IN) :: ulat(iim*jjm,llm)
345  REAL(rstd),INTENT(OUT)  :: un(3*iim*jjm,llm)
346
347  REAL(rstd)             :: uc(iim*jjm,llm,3)
348   
349    CALL compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc) 
350    CALL compute_wind_perp_from_wind_centered(uc, un)
351 
352 END SUBROUTINE compute_ulonlat2un
353
354
355END MODULE wind_mod
Note: See TracBrowser for help on using the repository browser.