source: codes/icosagcm/trunk/src/caldyn_sw.f90 @ 12

Last change on this file since 12 was 12, checked in by ymipsl, 12 years ago

dynamico tree creation

YM
File size: 20.1 KB
Line 
1MODULE caldyn_sw_mod
2  USE genmod
3  USE field_mod
4  USE transfert_mod
5  PRIVATE
6  TYPE(t_field),POINTER,SAVE :: f_Fe(:)
7  REAL(rstd),POINTER,SAVE    :: Fe(:)
8  TYPE(t_field),POINTER,SAVE :: f_Ki(:)
9  REAL(rstd),POINTER,SAVE    :: Ki(:)
10  TYPE(t_field),POINTER,SAVE :: f_qv(:)
11  REAL(rstd),POINTER,SAVE    :: qv(:)
12  TYPE(t_field),POINTER,SAVE :: f_qv2(:)
13  REAL(rstd),POINTER,SAVE    :: qv2(:)
14
15  TYPE(t_field),POINTER,SAVE :: f_t_tmp(:)
16  REAL(rstd),POINTER,SAVE    :: t_tmp(:)
17  TYPE(t_field),POINTER,SAVE :: f_u_tmp(:)
18  REAL(rstd),POINTER,SAVE    :: u_tmp(:)
19  TYPE(t_field),POINTER,SAVE :: f_z_tmp(:)
20  REAL(rstd),POINTER,SAVE    :: z_tmp(:)
21  TYPE(t_field),POINTER,SAVE :: f_dZ(:)
22  REAL(rstd),POINTER,SAVE    :: dZ(:)
23  TYPE(t_field),POINTER,SAVE :: f_diff_dhv(:)
24  REAL(rstd),POINTER,SAVE    :: diff_dhv(:)
25
26  TYPE(t_request),POINTER :: req(:) 
27  TYPE(t_request),POINTER :: req_u(:)   
28  PUBLIC :: allocate_caldyn,swap_caldyn,init_williamson,caldyn,write_caldyn,Compute_PV,Compute_enstrophy
29CONTAINS
30
31  SUBROUTINE allocate_caldyn
32  USE domain_mod
33  USE dimensions
34  USE geometry
35  USE metric
36  IMPLICIT NONE
37    INTEGER :: ind,i,j
38
39    CALL allocate_field(f_Fe,field_u,type_real)
40    CALL allocate_field(f_Ki,field_t,type_real)
41    CALL allocate_field(f_qv,field_z,type_real)
42    CALL allocate_field(f_qv2,field_z,type_real)
43    CALL allocate_field(f_t_tmp,field_t,type_real)
44    CALL allocate_field(f_u_tmp,field_u,type_real)
45    CALL allocate_field(f_z_tmp,field_z,type_real)
46    CALL allocate_field(f_dZ,field_z,type_real)
47    CALL allocate_field(f_diff_dhv,field_z,type_real)
48
49  CALL create_request(field_t,req)
50
51  DO ind=1,ndomain
52    DO i=ii_begin,ii_end+1
53      CALL request_add_point(ind,i,jj_begin-1,req)
54    ENDDO
55
56    DO j=jj_begin,jj_end
57      CALL request_add_point(ind,ii_end+1,j,req)
58    ENDDO   
59    DO i=ii_begin,ii_end
60      CALL request_add_point(ind,i,jj_end+1,req)
61    ENDDO   
62
63    DO j=jj_begin,jj_end+1
64      CALL request_add_point(ind,ii_begin-1,j,req)
65    ENDDO   
66   
67    DO i=ii_begin,ii_end
68      CALL request_add_point(ind,i,jj_begin,req)
69    ENDDO
70
71    DO j=jj_begin,jj_end
72      CALL request_add_point(ind,ii_end,j,req)
73    ENDDO   
74   
75    DO i=ii_begin,ii_end
76      CALL request_add_point(ind,i,jj_end,req)
77    ENDDO   
78
79    DO j=jj_begin,jj_end
80      CALL request_add_point(ind,ii_begin,j,req)
81    ENDDO   
82   
83  ENDDO
84
85  CALL create_request(field_u,req_u)
86  DO ind=1,ndomain
87    DO i=ii_begin,ii_end
88      CALL request_add_point(ind,i,jj_begin-1,req_u,rup)
89      CALL request_add_point(ind,i+1,jj_begin-1,req_u,lup)
90    ENDDO
91
92    DO j=jj_begin,jj_end
93      CALL request_add_point(ind,ii_end+1,j,req_u,left)
94      CALL request_add_point(ind,ii_end+1,j-1,req_u,lup)
95    ENDDO   
96   
97    DO i=ii_begin,ii_end
98      CALL request_add_point(ind,i,jj_end+1,req_u,ldown)
99      CALL request_add_point(ind,i-1,jj_end+1,req_u,rdown)
100    ENDDO   
101
102    DO j=jj_begin,jj_end
103      CALL request_add_point(ind,ii_begin-1,j,req_u,right)
104      CALL request_add_point(ind,ii_begin-1,j+1,req_u,rdown)
105    ENDDO   
106  ENDDO
107   
108  END SUBROUTINE allocate_caldyn
109 
110  SUBROUTINE swap_caldyn(ind)
111  IMPLICIT NONE
112    INTEGER,INTENT(IN) :: ind
113   
114      Fe=f_Fe(ind)
115      Ki=f_Ki(ind)
116      qv=f_qv(ind)
117      qv2=f_qv2(ind)
118      t_tmp=f_t_tmp(ind)
119      u_tmp=f_u_tmp(ind)
120      z_tmp=f_z_tmp(ind)
121      dZ=f_dZ(ind)
122      diff_dhv=f_diff_dhv(ind)
123     
124  END SUBROUTINE swap_caldyn
125
126 
127  SUBROUTINE caldyn(f_h, f_u, f_dh, f_du)
128  USE domain_mod
129  USE dimensions
130  USE grid_param
131  USE geometry
132  USE metric
133  USE write_field
134  IMPLICIT NONE
135  TYPE(t_field),POINTER :: f_h(:)
136  TYPE(t_field),POINTER :: f_u(:)
137  TYPE(t_field),POINTER :: f_dh(:)
138  TYPE(t_field),POINTER :: f_du(:)
139
140  REAL(rstd),POINTER :: h(:)
141  REAL(rstd),POINTER :: u(:)
142  REAL(rstd),POINTER :: dh(:)
143  REAL(rstd),POINTER :: du(:)
144  INTEGER :: ind
145  INTEGER,SAVE :: it=0
146 
147    CALL transfert_request(f_h,req_i1) 
148    CALL transfert_request(f_u,req_e1)
149    CALL transfert_request(f_u,req_e1) 
150   
151
152    DO ind=1,ndomain
153      CALL swap_dimensions(ind)
154      CALL swap_geometry(ind)
155      CALL swap_caldyn(ind)
156     
157      h=f_h(ind)
158      u=f_u(ind)
159      dh=f_dh(ind)
160      du=f_du(ind)
161     
162      CALL compute_caldyn(h, u, dh, du)
163
164    ENDDO
165
166
167
168    IF (mod(it,240)==0) THEN
169      CALL writefield("h",f_h)
170      CALL writefield("dh",f_dh)
171      CALL Compute_enstrophy
172    ENDIF
173    it=it+1     
174  END SUBROUTINE caldyn
175 
176 
177  SUBROUTINE compute_caldyn(hi,ue,dhi,due)
178  USE domain_mod
179  USE dimensions
180  USE geometry
181  USE metric
182  IMPLICIT NONE
183    REAL(rstd),INTENT(IN) :: hi(iim*jjm)
184    REAL(rstd),INTENT(IN) :: ue(iim*3*jjm)
185    REAL(rstd),INTENT(OUT):: dhi(iim*jjm)
186    REAL(rstd),INTENT(OUT):: due(iim*3*jjm)
187
188    REAL(rstd)            :: Fep(iim*3*jjm)
189   
190    INTEGER :: i,j,n
191    REAL(rstd) :: etav,hv
192   
193    DO j=jj_begin-1,jj_end+1
194      DO i=ii_begin-1,ii_end+1
195        n=(j-1)*iim+i
196       
197        Fe(n+u_right)=0.5*(hi(n)+hi(n+t_right))*ue(n+u_right)
198        Fe(n+u_lup)=0.5*(hi(n)+hi(n+t_lup))*ue(n+u_lup)
199        Fe(n+u_ldown)=0.5*(hi(n)+hi(n+t_ldown))*ue(n+u_ldown)
200     
201       u_tmp(n+u_right)=ue(n+u_right)
202       u_tmp(n+u_lup)=ue(n+u_lup)
203       u_tmp(n+u_ldown)=ue(n+u_ldown)
204       
205      ENDDO
206    ENDDO
207   
208   
209    DO j=jj_begin,jj_end
210      DO i=ii_begin,ii_end
211        n=(j-1)*iim+i
212
213        dhi(n)=-1./Ai(n)*(ne(n,right)*Fe(n+u_right)*le(n+u_right)  +  &
214                          ne(n,rup)*Fe(n+u_rup)*le(n+u_rup)        +  & 
215                          ne(n,lup)*Fe(n+u_lup)*le(n+u_lup)        +  & 
216                          ne(n,left)*Fe(n+u_left)*le(n+u_left)     +  & 
217                          ne(n,ldown)*Fe(n+u_ldown)*le(n+u_ldown)  +  & 
218                          ne(n,rdown)*Fe(n+u_rdown)*le(n+u_rdown))   
219
220      ENDDO
221    ENDDO
222
223    DO j=jj_begin,jj_end
224      DO i=ii_begin,ii_end
225        n=(j-1)*iim+i
226
227        Ki(n)=1/(4*Ai(n))*(le(n+u_right)*de(n+u_right)*ue(n+u_right)**2 +  &
228                           le(n+u_rup)*de(n+u_rup)*ue(n+u_rup)**2 +        &
229                           le(n+u_lup)*de(n+u_lup)*ue(n+u_lup)**2 +        &
230                           le(n+u_left)*de(n+u_left)*ue(n+u_left)**2 +     &
231                           le(n+u_ldown)*de(n+u_ldown)*ue(n+u_ldown)**2 +  &
232                           le(n+u_rdown)*de(n+u_rdown)*ue(n+u_rdown)**2 ) 
233       
234      ENDDO
235    ENDDO
236   
237    DO j=jj_begin,jj_end
238      DO i=ii_begin,ii_end
239        n=(j-1)*iim+i
240       
241        due(n+u_right)=1/de(n+u_right)*(ne(n,right)*(g*(hi(n)+bi(n))+Ki(n))+                            &
242                                        ne(n+t_right,left)*(g*(hi(n+t_right)+bi(n+t_right))+Ki(n+t_right)) )       
243       
244        due(n+u_lup)=1/de(n+u_lup)*(ne(n,lup)*(g*(hi(n)+bi(n))+Ki(n))+                                  &
245                                    ne(n+t_lup,rdown)*(g*(hi(n+t_lup)+bi(n+t_lup))+Ki(n+t_lup)) )       
246
247        due(n+u_ldown)=1/de(n+u_ldown)*(ne(n,ldown)*(g*(hi(n)+bi(n))+Ki(n))+                            &
248                                        ne(n+t_ldown,rup)*(g*(hi(n+t_ldown)+bi(n+t_ldown))+Ki(n+t_ldown)) )       
249      ENDDO
250    ENDDO
251
252   
253   
254    DO j=jj_begin-1,jj_end+1
255      DO i=ii_begin-1,ii_end+1
256        n=(j-1)*iim+i
257       
258          etav= 1./Av(n+z_up)*(  ne(n,rup)*ue(n+u_rup)*de(n+u_rup)                       &
259                               + ne(n+t_rup,left)*ue(n+t_rup+u_left)*de(n+t_rup+u_left)  &
260                               - ne(n,lup)*ue(n+u_lup)*de(n+u_lup) )                               
261
262          hv= Riv2(n,vup)*hi(n)+                  &
263              Riv2(n+t_rup,vldown)*hi(n+t_rup)+    &
264              Riv2(n+t_lup,vrdown)*hi(n+t_lup)
265     
266          qv(n+z_up)=(etav+fv(n+z_up))/hv
267          qv2(n+z_up)= (etav+fv(n+z_down))**2/hv
268          z_tmp(n+z_up)= Av(n+z_up)
269         
270          etav = 1./Av(n+z_down)*(  ne(n,ldown)*ue(n+u_ldown)*de(n+u_ldown)                          &
271                                  + ne(n+t_ldown,right)*ue(n+t_ldown+u_right)*de(n+t_ldown+u_right)  &
272                                  - ne(n,rdown)*ue(n+u_rdown)*de(n+u_rdown) )
273       
274          hv = Riv2(n,vdown)*hi(n)+                  &
275               Riv2(n+t_ldown,vrup)*hi(n+t_ldown)+   &
276               Riv2(n+t_rdown,vlup)*hi(n+t_rdown)
277                       
278          qv(n+z_down)=(etav+fv(n+z_down))/hv
279          qv2(n+z_down)= (etav+fv(n+z_down))**2/hv
280         
281          z_tmp(n+z_down)=Av(n+z_down)
282        ENDDO
283      ENDDO
284
285    DO j=jj_begin,jj_end
286      DO i=ii_begin,ii_end
287        n=(j-1)*iim+i
288
289        due(n+u_right) = due(n+u_right)+0.5*(qv(n+z_rdown)+qv(n+z_rup))/de(n+u_right) *     & 
290                       ( wee(n+u_right,1,1)*le(n+u_rup)*Fe(n+u_rup)+                        &
291                         wee(n+u_right,2,1)*le(n+u_lup)*Fe(n+u_lup)+                        &
292                         wee(n+u_right,3,1)*le(n+u_left)*Fe(n+u_left)+                      &
293                         wee(n+u_right,4,1)*le(n+u_ldown)*Fe(n+u_ldown)+                    &
294                         wee(n+u_right,5,1)*le(n+u_rdown)*Fe(n+u_rdown)+                    & 
295                         wee(n+u_right,1,2)*le(n+t_right+u_ldown)*Fe(n+t_right+u_ldown)+    &
296                         wee(n+u_right,2,2)*le(n+t_right+u_rdown)*Fe(n+t_right+u_rdown)+    &
297                         wee(n+u_right,3,2)*le(n+t_right+u_right)*Fe(n+t_right+u_right)+    &
298                         wee(n+u_right,4,2)*le(n+t_right+u_rup)*Fe(n+t_right+u_rup)+        &
299                         wee(n+u_right,5,2)*le(n+t_right+u_lup)*Fe(n+t_right+u_lup) )
300 
301
302        due(n+u_lup) = due(n+u_lup) + 0.5*(qv(n+z_up)+qv(n+z_lup))/de(n+u_lup) *            & 
303                       ( wee(n+u_lup,1,1)*le(n+u_left)*Fe(n+u_left)+                        &
304                         wee(n+u_lup,2,1)*le(n+u_ldown)*Fe(n+u_ldown)+                      &
305                         wee(n+u_lup,3,1)*le(n+u_rdown)*Fe(n+u_rdown)+                      &
306                         wee(n+u_lup,4,1)*le(n+u_right)*Fe(n+u_right)+                      &
307                         wee(n+u_lup,5,1)*le(n+u_rup)*Fe(n+u_rup)+                          & 
308                         wee(n+u_lup,1,2)*le(n+t_lup+u_right)*Fe(n+t_lup+u_right)+          &
309                         wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*Fe(n+t_lup+u_rup)+              &
310                         wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*Fe(n+t_lup+u_lup)+              &
311                         wee(n+u_lup,4,2)*le(n+t_lup+u_left)*Fe(n+t_lup+u_left)+            &
312                         wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*Fe(n+t_lup+u_ldown) )
313
314
315        due(n+u_ldown) = due(n+u_ldown)+ 0.5*(qv(n+z_ldown)+qv(n+z_down))/de(n+u_ldown) *   & 
316                       ( wee(n+u_ldown,1,1)*le(n+u_rdown)*Fe(n+u_rdown)+                    &
317                         wee(n+u_ldown,2,1)*le(n+u_right)*Fe(n+u_right)+                    &
318                         wee(n+u_ldown,3,1)*le(n+u_rup)*Fe(n+u_rup)+                        &
319                         wee(n+u_ldown,4,1)*le(n+u_lup)*Fe(n+u_lup)+                        &
320                         wee(n+u_ldown,5,1)*le(n+u_left)*Fe(n+u_left)+                      & 
321                         wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*Fe(n+t_ldown+u_lup)+        &
322                         wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*Fe(n+t_ldown+u_left)+      &
323                         wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*Fe(n+t_ldown+u_ldown)+    &
324                         wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*Fe(n+t_ldown+u_rdown)+    &
325                         wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*Fe(n+t_ldown+u_right) )
326
327     
328      ENDDO
329    ENDDO
330   
331    DO j=jj_begin-1,jj_end+1
332      DO i=ii_begin-1,ii_end+1
333        n=(j-1)*iim+i
334       
335          dZ(n+z_up) = 2*qv(n+z_up)/Av(n+z_up)*(  ne(n,rup)*due(n+u_rup)*de(n+u_rup)                       &
336                                                + ne(n+t_rup,left)*due(n+t_rup+u_left)*de(n+t_rup+u_left)  &
337                                                - ne(n,lup)*due(n+u_lup)*de(n+u_lup) )
338          dZ(n+z_up)= dZ(n+z_up)-(Riv2(n,vup)*dhi(n)+                  &
339                                   Riv2(n+t_rup,vldown)*dhi(n+t_rup)+    &
340                                   Riv2(n+t_lup,vrdown)*dhi(n+t_lup))*qv(n+z_up)**2
341                                             
342          dZ(n+z_down) = 2*qv(n+z_down)/Av(n+z_down)*(  ne(n,ldown)*due(n+u_ldown)*de(n+u_ldown)                          &
343                                                      + ne(n+t_ldown,right)*due(n+t_ldown+u_right)*de(n+t_ldown+u_right)  &
344                                                      - ne(n,rdown)*due(n+u_rdown)*de(n+u_rdown) )
345          dZ(n+z_down)= dZ(n+z_down)- (Riv2(n,vdown)*dhi(n)+                  &
346                                      Riv2(n+t_ldown,vrup)*dhi(n+t_ldown)+   &
347                                      Riv2(n+t_rdown,vlup)*dhi(n+t_rdown))*qv(n+z_down)**2
348                                                     
349      ENDDO
350     ENDDO
351
352    DO j=jj_begin,jj_end
353      DO i=ii_begin,ii_end
354        n=(j-1)*iim+i
355
356        Fep(n+u_right) = 1./de(n+u_right) *     & 
357                       ( wee(n+u_right,1,1)*le(n+u_rup)*Fe(n+u_rup)+                        &
358                         wee(n+u_right,2,1)*le(n+u_lup)*Fe(n+u_lup)+                        &
359                         wee(n+u_right,3,1)*le(n+u_left)*Fe(n+u_left)+                      &
360                         wee(n+u_right,4,1)*le(n+u_ldown)*Fe(n+u_ldown)+                    &
361                         wee(n+u_right,5,1)*le(n+u_rdown)*Fe(n+u_rdown)+                    & 
362                         wee(n+u_right,1,2)*le(n+t_right+u_ldown)*Fe(n+t_right+u_ldown)+    &
363                         wee(n+u_right,2,2)*le(n+t_right+u_rdown)*Fe(n+t_right+u_rdown)+    &
364                         wee(n+u_right,3,2)*le(n+t_right+u_right)*Fe(n+t_right+u_right)+    &
365                         wee(n+u_right,4,2)*le(n+t_right+u_rup)*Fe(n+t_right+u_rup)+        &
366                         wee(n+u_right,5,2)*le(n+t_right+u_lup)*Fe(n+t_right+u_lup) )
367 
368
369        Fep(n+u_lup) = 1/de(n+u_lup) *            & 
370                       ( wee(n+u_lup,1,1)*le(n+u_left)*Fe(n+u_left)+                        &
371                         wee(n+u_lup,2,1)*le(n+u_ldown)*Fe(n+u_ldown)+                      &
372                         wee(n+u_lup,3,1)*le(n+u_rdown)*Fe(n+u_rdown)+                      &
373                         wee(n+u_lup,4,1)*le(n+u_right)*Fe(n+u_right)+                      &
374                         wee(n+u_lup,5,1)*le(n+u_rup)*Fe(n+u_rup)+                          & 
375                         wee(n+u_lup,1,2)*le(n+t_lup+u_right)*Fe(n+t_lup+u_right)+          &
376                         wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*Fe(n+t_lup+u_rup)+              &
377                         wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*Fe(n+t_lup+u_lup)+              &
378                         wee(n+u_lup,4,2)*le(n+t_lup+u_left)*Fe(n+t_lup+u_left)+            &
379                         wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*Fe(n+t_lup+u_ldown) )
380
381
382        Fep(n+u_ldown) = 1./de(n+u_ldown) *   & 
383                       ( wee(n+u_ldown,1,1)*le(n+u_rdown)*Fe(n+u_rdown)+                    &
384                         wee(n+u_ldown,2,1)*le(n+u_right)*Fe(n+u_right)+                    &
385                         wee(n+u_ldown,3,1)*le(n+u_rup)*Fe(n+u_rup)+                        &
386                         wee(n+u_ldown,4,1)*le(n+u_lup)*Fe(n+u_lup)+                        &
387                         wee(n+u_ldown,5,1)*le(n+u_left)*Fe(n+u_left)+                      & 
388                         wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*Fe(n+t_ldown+u_lup)+        &
389                         wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*Fe(n+t_ldown+u_left)+      &
390                         wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*Fe(n+t_ldown+u_ldown)+    &
391                         wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*Fe(n+t_ldown+u_rdown)+    &
392                         wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*Fe(n+t_ldown+u_right) )
393
394     
395      ENDDO
396    ENDDO
397
398    DO j=jj_begin-1,jj_end+1
399      DO i=ii_begin-1,ii_end+1
400        n=(j-1)*iim+i
401         
402          diff_dhv(n+z_up)= 1./Av(n+z_up)*(  ne(n,rup)*Fep(n+u_rup)*de(n+u_rup)                       &
403                                           + ne(n+t_rup,left)*Fep(n+t_rup+u_left)*de(n+t_rup+u_left)  &
404                                           - ne(n,lup)*Fep(n+u_lup)*de(n+u_lup) )                               
405
406          diff_dhv(n+z_up)=diff_dhv(n+z_up)-(Riv2(n,vup)*dhi(n)+                   &
407                                             Riv2(n+t_rup,vldown)*dhi(n+t_rup)+    &
408                                             Riv2(n+t_lup,vrdown)*dhi(n+t_lup))
409         
410         
411          diff_dhv(n+z_down) = 1./Av(n+z_down)*(  ne(n,ldown)*Fep(n+u_ldown)*de(n+u_ldown)                          &
412                                                + ne(n+t_ldown,right)*Fep(n+t_ldown+u_right)*de(n+t_ldown+u_right)  &
413                                                - ne(n,rdown)*Fep(n+u_rdown)*de(n+u_rdown) )
414
415          diff_dhv(n+z_down) = diff_dhv(n+z_down) - (Riv2(n,vdown)*dhi(n)+                  &
416                                                     Riv2(n+t_ldown,vrup)*dhi(n+t_ldown)+   &
417                                                     Riv2(n+t_rdown,vlup)*dhi(n+t_rdown))
418       ENDDO
419    ENDDO   
420                                                                         
421  END SUBROUTINE compute_caldyn
422   
423  SUBROUTINE write_caldyn
424  USE write_field
425  USE domain_mod
426  USE dimensions
427  USE geometry
428  USE metric
429  IMPLICIT NONE
430   
431    INTEGER :: ind,i,j,n
432   
433    CALL transfert_request(f_u_tmp,req_u)
434
435    DO ind=1,ndomain
436      CALL swap_caldyn(ind)
437      CALL swap_geometry(ind)
438      CALL swap_dimensions(ind)
439     
440      DO j=jj_begin,jj_end
441        DO i=ii_begin,ii_end
442          n=(j-1)*iim+i
443
444          Ki(n)=1/(4*Ai(n))*(le(n+u_right)*de(n+u_right)*u_tmp(n+u_right)**2 +  &
445                             le(n+u_rup)*de(n+u_rup)*u_tmp(n+u_rup)**2 +        &
446                             le(n+u_lup)*de(n+u_lup)*u_tmp(n+u_lup)**2 +        &
447                             le(n+u_left)*de(n+u_left)*u_tmp(n+u_left)**2 +     &
448                             le(n+u_ldown)*de(n+u_ldown)*u_tmp(n+u_ldown)**2 +  &
449                             le(n+u_rdown)*de(n+u_rdown)*u_tmp(n+u_rdown)**2 ) 
450       
451        ENDDO
452      ENDDO
453    ENDDO
454   
455    CALL writefield("Ki",f_Ki)
456
457  END SUBROUTINE write_caldyn
458     
459  SUBROUTINE Compute_PV
460  USE domain_mod
461  USE dimensions
462  USE geometry
463  USE metric
464  USE write_field
465  IMPLICIT NONE
466    REAL(rstd) :: PV
467    INTEGER :: i,j,n,ind
468    PV=0
469    DO ind=1,ndomain
470      CALL swap_caldyn(ind)
471      CALL swap_geometry(ind)
472      CALL swap_dimensions(ind)
473 
474      DO j=jj_begin+1,jj_end
475        DO i=ii_begin,ii_end-1
476          n=(j-1)*iim+i
477          qv(n+z_down)=qv(n+z_down)*Av(n+z_down)
478          PV=PV+qv(n+z_down)
479        ENDDO
480      ENDDO
481
482      DO j=jj_begin,jj_end-1
483        DO i=ii_begin+1,ii_end
484          n=(j-1)*iim+i
485          qv(n+z_up)=qv(n+z_up)*Av(n+z_up)
486          PV=PV+qv(n+z_up)
487        ENDDO
488      ENDDO
489    ENDDO
490    CALL writefield("PV",f_qv)
491   
492    PRINT *,"PV=",PV 
493  END SUBROUTINE Compute_PV   
494 
495  SUBROUTINE Compute_enstrophy
496  USE domain_mod
497  USE dimensions
498  USE geometry
499  USE metric
500  USE write_field
501 
502  IMPLICIT NONE
503    REAL(rstd) :: Es
504    INTEGER :: i,j,n,ind
505    Es=0
506    DO ind=1,ndomain
507      CALL swap_caldyn(ind)
508      CALL swap_geometry(ind)
509      CALL swap_dimensions(ind)
510 
511      DO j=jj_begin+1,jj_end
512        DO i=ii_begin,ii_end-1
513          n=(j-1)*iim+i
514          Es=Es+qv2(n+z_down)*Av(n+z_down)
515        ENDDO
516      ENDDO
517
518      DO j=jj_begin,jj_end-1
519        DO i=ii_begin+1,ii_end
520          n=(j-1)*iim+i
521          Es=Es+qv2(n+z_up)*Av(n+z_up)
522        ENDDO
523      ENDDO
524    ENDDO
525   
526    PRINT *,"Enstrophy=",Es 
527    CALL writefield("Ens",f_qv2)
528
529    Es=0
530    DO ind=1,ndomain
531      CALL swap_caldyn(ind)
532      CALL swap_geometry(ind)
533      CALL swap_dimensions(ind)
534 
535      DO j=jj_begin+1,jj_end
536        DO i=ii_begin,ii_end-1
537          n=(j-1)*iim+i
538          Es=Es+dZ(n+z_down)*Av(n+z_down)
539        ENDDO
540      ENDDO
541
542      DO j=jj_begin,jj_end-1
543        DO i=ii_begin+1,ii_end
544          n=(j-1)*iim+i
545          Es=Es+dZ(n+z_up)*Av(n+z_up)
546        ENDDO
547      ENDDO
548    ENDDO
549   
550    PRINT *,"D/dt Enstrophy=",Es 
551    CALL writefield("dEns",f_dZ)
552    CALL writefield("diff_dhv",f_diff_dhv)
553         
554  END SUBROUTINE Compute_enstrophy
555 
556END MODULE caldyn_sw_mod
Note: See TracBrowser for help on using the repository browser.