source: codes/icosagcm/trunk/src/dissip_sw.f90 @ 176

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