source: codes/icosagcm/trunk/src/dissip.f90 @ 13

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

some update

YM
File size: 11.6 KB
Line 
1MODULE dissip_mod
2  USE genmod
3  USE field_mod
4  USE transfert_mod
5
6  TYPE(t_field),POINTER,SAVE :: f_gradrot(:)
7  TYPE(t_request),POINTER :: req_dissip(:) 
8  TYPE(t_field),POINTER,SAVE :: f_due(:)
9 
10  INTEGER,PARAMETER :: nitergdiv=1
11  INTEGER,PARAMETER :: nitergrot=1
12
13  REAL :: tetagdiv
14  REAL :: tetagrot
15 
16  REAL :: cdivu
17  REAL :: crot
18  INTEGER :: idissip
19  REAL    :: dtdissip
20 
21CONTAINS
22
23  SUBROUTINE allocate_dissip
24  IMPLICIT NONE 
25    CALL allocate_field(f_gradrot,field_u,type_real)
26    CALL allocate_field(f_due,field_u,type_real)
27  END SUBROUTINE allocate_dissip
28 
29  SUBROUTINE init_dissip(dt)
30  USE domain_mod
31  USE dimensions
32  USE geometry
33  USE metric
34  USE ioipsl
35 
36  IMPLICIT NONE
37  REAL,INTENT(IN)       :: dt
38 
39  TYPE(t_field),POINTER :: f_u(:)
40  TYPE(t_field),POINTER :: f_du(:)
41  REAL(rstd),POINTER    :: u(:)
42  REAL(rstd),POINTER    :: du(:)
43  REAL(rstd)            :: dumax,dumaxmax,dumaxm1
44  REAL(rstd)            :: r
45  INTEGER :: i,j,n,ind,it
46  REAL :: sum1,sum2
47
48    tetagdiv=5000
49    CALL getin("tetagdiv",tetagdiv)
50    tetagrot=5000
51    CALL getin("tetagrot",tetagrot)
52   
53    CALL allocate_dissip
54 
55    CALL allocate_field(f_u,field_u,type_real)
56    CALL allocate_field(f_du,field_u,type_real)
57    CALL create_request(field_u,req_dissip)
58
59    DO ind=1,ndomain
60      DO i=ii_begin,ii_end
61        CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup)
62        CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup)
63      ENDDO
64
65      DO j=jj_begin,jj_end
66        CALL request_add_point(ind,ii_end+1,j,req_dissip,left)
67        CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup)
68      ENDDO   
69   
70      DO i=ii_begin,ii_end
71        CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown)
72        CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown)
73      ENDDO   
74
75      DO j=jj_begin,jj_end
76        CALL request_add_point(ind,ii_begin-1,j,req_dissip,right)
77        CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown)
78      ENDDO   
79
80      DO i=ii_begin+1,ii_end-1
81        CALL request_add_point(ind,i,jj_begin,req_dissip,right)
82        CALL request_add_point(ind,i,jj_end,req_dissip,right)
83      ENDDO
84   
85      DO j=jj_begin+1,jj_end-1
86        CALL request_add_point(ind,ii_begin,j,req_dissip,rup)
87        CALL request_add_point(ind,ii_end,j,req_dissip,rup)
88      ENDDO   
89
90      CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left)
91      CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown)
92      CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left)
93      CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown)
94     
95    ENDDO
96
97
98    cdivu=-1
99    crot=-1
100 
101    CALL RANDOM_SEED()
102   
103    DO ind=1,ndomain
104      CALL swap_dimensions(ind)
105      CALL swap_geometry(ind)
106      u=f_u(ind)
107
108      DO j=jj_begin,jj_end
109        DO i=ii_begin,ii_end
110          n=(j-1)*iim+i
111          CALL RANDOM_NUMBER(r)
112          u(n+u_right)=r-0.5
113          CALL RANDOM_NUMBER(r)
114          u(n+u_lup)=r-0.5
115          CALL RANDOM_NUMBER(r)
116          u(n+u_ldown)=r-0.5
117       ENDDO
118      ENDDO
119       
120    ENDDO
121 
122
123    CALL transfert_request(f_u,req_dissip)
124    CALL transfert_request(f_u,req_dissip)
125
126    DO it=1,500
127
128      dumaxm1=dumax 
129      dumax=0
130      DO ind=1,ndomain
131        CALL swap_dimensions(ind)
132        CALL swap_geometry(ind)
133        u=f_u(ind)
134        du=f_du(ind)
135        CALL gradiv(u,du)
136      ENDDO
137      CALL transfert_request(f_du,req_dissip)
138      CALL transfert_request(f_du,req_dissip)
139     
140      DO ind=1,ndomain
141        CALL swap_dimensions(ind)
142        CALL swap_geometry(ind)
143        du=f_du(ind)
144       
145        DO j=jj_begin,jj_end
146          DO i=ii_begin,ii_end
147            n=(j-1)*iim+i
148            if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
149            if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
150            if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
151          ENDDO
152        ENDDO
153      ENDDO
154      u=du/dumax
155      PRINT *,"gradiv : it :",it ,": dumax",(dumax+dumaxm1)/2
156
157    ENDDO 
158    PRINT *,"gradiv : dumax",(dumax+dumaxm1)/2
159   
160    cdivu=dumax**(-1./nitergdiv)
161    PRINT *,"cdivu : ",cdivu
162
163    DO ind=1,ndomain
164      CALL swap_dimensions(ind)
165      CALL swap_geometry(ind)
166      u=f_u(ind)
167
168     DO j=jj_begin,jj_end
169        DO i=ii_begin,ii_end
170          n=(j-1)*iim+i
171          CALL RANDOM_NUMBER(r)
172          u(n+u_right)=r-0.5
173          CALL RANDOM_NUMBER(r)
174          u(n+u_lup)=r-0.5
175          CALL RANDOM_NUMBER(r)
176          u(n+u_ldown)=r-0.5
177       ENDDO
178      ENDDO
179       
180    ENDDO
181
182    CALL transfert_request(f_u,req_dissip)
183    CALL transfert_request(f_u,req_dissip)
184 
185    DO it=1,500
186
187      dumaxm1=dumax 
188      dumax=0
189      sum1=0 ; sum2=0
190      DO ind=1,ndomain
191        u=f_u(ind)
192        du=f_du(ind)
193        CALL gradrot(u,du,ind,sum1,sum2)
194      ENDDO
195
196      CALL transfert_request(f_du,req_dissip)
197      CALL transfert_request(f_du,req_dissip)
198     
199      DO ind=1,ndomain
200        CALL swap_dimensions(ind)
201        CALL swap_geometry(ind)
202        du=f_du(ind)
203       
204        DO j=jj_begin,jj_end
205          DO i=ii_begin,ii_end
206            n=(j-1)*iim+i
207            if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
208            if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
209            if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
210          ENDDO
211        ENDDO
212      ENDDO
213      u=du/dumax
214      PRINT *,"gradrot : it :",it ,": dumax",(dumax+dumaxm1)/2
215
216    ENDDO 
217    PRINT *,"gradrot : dumax",(dumax+dumaxm1)/2
218 
219    crot=dumax**(-1/nitergrot) 
220    PRINT *,"crot : ",crot
221   
222    idissip=INT(MIN(tetagdiv,tetagrot)/(2.*dt))
223    idissip=MAX(1,idissip)
224    dtdissip=idissip*dt
225    PRINT *,"idissip",idissip," dtdissip ",dtdissip
226
227  END SUBROUTINE init_dissip 
228 
229 
230  SUBROUTINE dissip(f_ue)
231  USE domain_mod
232  USE dimensions
233  USE geometry
234  USE metric
235  IMPLICIT NONE
236    TYPE(t_field),POINTER :: f_ue(:)
237!    TYPE(t_field),POINTER :: f_due(:)
238    REAL(rstd),POINTER         :: ue(:)
239    REAL(rstd),POINTER         :: due(:)
240    REAL(rstd),POINTER         :: gradrot_e(:)
241    INTEGER :: ind
242    REAL :: v
243    REAL :: sum1,sum2
244   
245    CALL transfert_request(f_ue,req_dissip)
246    CALL transfert_request(f_ue,req_dissip)
247    sum1=0 ; sum2=0
248    DO ind=1,ndomain
249      CALL swap_dimensions(ind)
250      CALL swap_geometry(ind)
251      ue=f_ue(ind)
252      due=f_due(ind) 
253      gradrot_e=f_gradrot(ind)
254      CALL gradiv(ue,due)
255      due=due*dtdissip/tetagdiv
256      CALL gradrot(ue,gradrot_e,ind,sum1,sum2)
257      due=due+gradrot_e*dtdissip/tetagrot 
258    ENDDO
259    PRINT *,"SUM1 == SUM2 ??  ", sum1,sum2
260    DO ind=1,ndomain
261      CALL swap_dimensions(ind)
262      CALL swap_geometry(ind)
263      ue=f_ue(ind)
264      due=f_due(ind) 
265     
266      ue(:)=ue(:)+0.5*due(:)
267    ENDDO
268     
269  END SUBROUTINE dissip
270 
271  SUBROUTINE gradiv(ue,gradivu_e)
272  USE domain_mod
273  USE dimensions
274  USE geometry
275  USE metric
276  IMPLICIT NONE
277    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm)
278    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm)
279    REAL(rstd)             :: divu_i(iim*jjm)
280
281    INTEGER :: i,j,n
282   
283    DO j=jj_begin,jj_end
284      DO i=ii_begin,ii_end
285        n=(j-1)*iim+i
286        divu_i(n)=1./Ai(n)*(ne(n,right)*ue(n+u_right)*le(n+u_right)  +  &
287                             ne(n,rup)*ue(n+u_rup)*le(n+u_rup)        +  & 
288                             ne(n,lup)*ue(n+u_lup)*le(n+u_lup)        +  & 
289                             ne(n,left)*ue(n+u_left)*le(n+u_left)     +  & 
290                             ne(n,ldown)*ue(n+u_ldown)*le(n+u_ldown)  +  & 
291                             ne(n,rdown)*ue(n+u_rdown)*le(n+u_rdown))
292      ENDDO
293    ENDDO
294   
295    DO j=jj_begin,jj_end
296      DO i=ii_begin,ii_end
297 
298        n=(j-1)*iim+i
299 
300        gradivu_e(n+u_right)=-1/de(n+u_right)*(ne(n,right)*divu_i(n)+ ne(n+t_right,left)*divu_i(n+t_right) )       
301
302        gradivu_e(n+u_lup)=-1/de(n+u_lup)*(ne(n,lup)*divu_i(n)+ ne(n+t_lup,rdown)*divu_i(n+t_lup ))       
303   
304        gradivu_e(n+u_ldown)=-1/de(n+u_ldown)*(ne(n,ldown)*divu_i(n)+ne(n+t_ldown,rup)*divu_i(n+t_ldown) )
305
306      ENDDO
307    ENDDO
308
309    DO j=jj_begin,jj_end
310      DO i=ii_begin,ii_end
311        n=(j-1)*iim+i
312        gradivu_e(n+u_right)=gradivu_e(n+u_right)*cdivu
313        gradivu_e(n+u_lup)=gradivu_e(n+u_lup)*cdivu
314        gradivu_e(n+u_ldown)=gradivu_e(n+u_ldown)*cdivu
315      ENDDO
316    ENDDO
317
318
319   
320  END SUBROUTINE gradiv
321 
322 
323  SUBROUTINE gradrot(ue,gradrot_e,ind,sum1,sum2)
324  USE domain_mod
325  USE dimensions
326  USE geometry
327  USE metric
328  IMPLICIT NONE
329    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm)
330    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm)
331    INTEGER,INTENT(IN) :: ind
332    REAL,INTENT(OUT) ::sum1,sum2
333    REAL(rstd)             :: rot_v(2*iim*jjm)
334
335    INTEGER :: i,j,n
336 
337    DO j=jj_begin-1,jj_end+1
338      DO i=ii_begin-1,ii_end+1
339        n=(j-1)*iim+i
340       
341        rot_v(n+z_up)= 1./Av(n+z_up)*(  ne(n,rup)*ue(n+u_rup)*de(n+u_rup)                       &
342                              + ne(n+t_rup,left)*ue(n+t_rup+u_left)*de(n+t_rup+u_left)          &
343                              - ne(n,lup)*ue(n+u_lup)*de(n+u_lup) ) 
344                             
345        rot_v(n+z_down) = 1./Av(n+z_down)*(   ne(n,ldown)*ue(n+u_ldown)*de(n+u_ldown)                &
346                                  + ne(n+t_ldown,right)*ue(n+t_ldown+u_right)*de(n+t_ldown+u_right)  &
347                                  - ne(n,rdown)*ue(n+u_rdown)*de(n+u_rdown) )
348
349      ENDDO
350    ENDDO                             
351 
352   
353    DO j=jj_begin,jj_end
354      DO i=ii_begin,ii_end
355        n=(j-1)*iim+i
356 
357        gradrot_e(n+u_right)=1/le(n+u_right)*ne(n,right)*(rot_v(n+z_rdown)-rot_v(n+z_rup)) 
358       
359        gradrot_e(n+u_lup)=1/le(n+u_lup)*ne(n,lup)*(rot_v(n+z_up)-rot_v(n+z_lup)) 
360       
361        gradrot_e(n+u_ldown)=1/le(n+u_ldown)*ne(n,ldown)*(rot_v(n+z_ldown)-rot_v(n+z_down))
362       
363      ENDDO
364    ENDDO
365
366   
367    DO j=jj_begin,jj_end-1
368      DO i=ii_begin,ii_end-1
369        n=(j-1)*iim+i
370        sum1=sum1+rot_v(n+z_rup)**2*Av(n+z_rup)
371      ENDDO
372    ENDDO
373   
374    DO j=jj_begin+1,jj_end
375      DO i=ii_begin,ii_end-1
376        n=(j-1)*iim+i
377        sum1=sum1+rot_v(n+z_rdown)**2*Av(n+z_rdown)
378      ENDDO
379    ENDDO
380   
381   
382    DO j=jj_begin+1,jj_end-1
383      DO i=ii_begin,ii_end-1
384        n=(j-1)*iim+i
385        sum2=sum2+gradrot_e(n+u_right)*le(n+u_right)*de(n+u_right)*ue(n+u_right)
386      ENDDO
387    ENDDO
388   
389    DO j=jj_begin,jj_end-1
390      DO i=ii_begin+1,ii_end-1
391        n=(j-1)*iim+i
392        sum2=sum2+gradrot_e(n+u_rup)*le(n+u_rup)*de(n+u_rup)*ue(n+u_rup)
393      ENDDO
394    ENDDO
395
396    DO j=jj_begin,jj_end-1
397      DO i=ii_begin+1,ii_end
398        n=(j-1)*iim+i
399        sum2=sum2+gradrot_e(n+u_lup)*le(n+u_lup)*de(n+u_lup)*ue(n+u_lup)
400      ENDDO
401    ENDDO
402
403! frontière
404      j=jj_begin
405      DO i=ii_begin,ii_end-1
406        n=(j-1)*iim+i
407        if (domain(ind)%own(i,j)) sum2=sum2+gradrot_e(n+u_right)*le(n+u_right)*de(n+u_right)*ue(n+u_right)
408      ENDDO
409     
410      j=jj_end
411      DO i=ii_begin,ii_end-1
412        n=(j-1)*iim+i
413        if (domain(ind)%own(i,j)) sum2=sum2+gradrot_e(n+u_right)*le(n+u_right)*de(n+u_right)*ue(n+u_right)
414      ENDDO
415     
416     i=ii_begin
417     DO j=jj_begin,jj_end-1
418        n=(j-1)*iim+i
419        if (domain(ind)%own(i,j))sum2=sum2+gradrot_e(n+u_rup)*le(n+u_rup)*de(n+u_rup)*ue(n+u_rup)
420     ENDDO
421
422     i=ii_end
423     DO j=jj_begin,jj_end-1
424        n=(j-1)*iim+i
425        if (domain(ind)%own(i,j))sum2=sum2+gradrot_e(n+u_rup)*le(n+u_rup)*de(n+u_rup)*ue(n+u_rup)
426     ENDDO
427
428     
429     
430     
431    DO j=jj_begin,jj_end
432      DO i=ii_begin,ii_end
433        n=(j-1)*iim+i
434 
435        gradrot_e(n+u_right)=-gradrot_e(n+u_right)*crot       
436        gradrot_e(n+u_lup)=-gradrot_e(n+u_lup)*crot
437        gradrot_e(n+u_ldown)=-gradrot_e(n+u_ldown)*crot
438       
439      ENDDO
440    ENDDO 
441   
442   
443  END SUBROUTINE
444
445
446END MODULE dissip_mod
447       
Note: See TracBrowser for help on using the repository browser.