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

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

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