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

Last change on this file since 155 was 146, checked in by ymipsl, 11 years ago

Set constant sign for wind way :
ne(ij,right)==ne_right=1
ne(ij,rup)==ne_rup=-1
ne(ij,lup)==ne_lup=1
ne(ij,left)==ne_left=-1
ne(ij,ldown)==ne_ldown=1
ne(ij,rdown)==ne_rdown=-1

Modified transfert function to be compliant for this convention.

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