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

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

dynamico tree creation

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