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

Last change on this file since 67 was 44, checked in by dubos, 12 years ago

Small fixes to enable compilation on MacOSX
Introduced key CPP_NETCDF4 (netcdf_mod.F90) to enable NetCDF4 with sequential version.

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