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

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

Simplify the management of the module.

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