source: codes/icosagcm/trunk/src/dissip_gcm.f90 @ 22

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

Review of 3D transport code
Fewer interations in dissip_gcm
This revision is certainly broken ; see next revision
File size: 17.9 KB
Line 
1MODULE dissip_gcm_mod
2  USE icosa
3
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  INTEGER,PARAMETER :: niterdivgrad=1
11
12  REAL,ALLOCATABLE,SAVE :: tau_graddiv(:)
13  REAL,ALLOCATABLE,SAVE :: tau_gradrot(:)
14  REAL,ALLOCATABLE,SAVE :: tau_divgrad(:)
15 
16  REAL,SAVE :: cgraddiv
17  REAL,SAVE :: cgradrot
18  REAL,SAVE :: cdivgrad
19 
20  INTEGER :: idissip
21  REAL    :: dtdissip
22 
23CONTAINS
24
25  SUBROUTINE allocate_dissip
26  USE icosa
27  IMPLICIT NONE 
28    CALL allocate_field(f_gradrot,field_u,type_real,llm)
29    CALL allocate_field(f_due,field_u,type_real,llm)
30    ALLOCATE(tau_graddiv(llm))
31    ALLOCATE(tau_gradrot(llm))   
32    ALLOCATE(tau_divgrad(llm))
33  END SUBROUTINE allocate_dissip
34 
35  SUBROUTINE init_dissip(dt)
36  USE icosa
37  USE disvert_mod
38 
39  IMPLICIT NONE
40  REAL,INTENT(IN)       :: dt
41 
42  TYPE(t_field),POINTER :: f_u(:)
43  TYPE(t_field),POINTER :: f_du(:)
44  REAL(rstd),POINTER    :: u(:)
45  REAL(rstd),POINTER    :: du(:)
46  TYPE(t_field),POINTER :: f_theta(:)
47  TYPE(t_field),POINTER :: f_dtheta(:)
48  REAL(rstd),POINTER    :: theta(:)
49  REAL(rstd),POINTER    :: dtheta(:)
50  REAL(rstd)            :: dumax
51  REAL(rstd)            :: dthetamax
52  REAL(rstd)            :: r
53  REAL(rstd)            :: tau
54  REAL(rstd)            :: zz, zvert, fact
55  INTEGER               :: l
56 
57           
58  INTEGER :: i,j,n,ind,it
59
60    CALL allocate_dissip
61    CALL allocate_field(f_u,field_u,type_real)
62    CALL allocate_field(f_du,field_u,type_real)
63    CALL allocate_field(f_theta,field_t,type_real)
64    CALL allocate_field(f_dtheta,field_t,type_real)
65
66    tau_graddiv(:)=5000
67    CALL getin("tau_graddiv",tau)
68    tau_graddiv(:)=tau
69   
70    tau_gradrot(:)=5000
71    CALL getin("tau_gradrot",tau)
72    tau_gradrot(:)=tau
73   
74
75    tau_divgrad(:)=5000
76    CALL getin("tau_divgrad",tau)
77    tau_divgrad(:)=tau
78 
79    CALL create_request(field_u,req_dissip)
80
81    DO ind=1,ndomain
82      DO i=ii_begin,ii_end
83        CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup)
84        CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup)
85      ENDDO
86
87      DO j=jj_begin,jj_end
88        CALL request_add_point(ind,ii_end+1,j,req_dissip,left)
89        CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup)
90      ENDDO   
91   
92      DO i=ii_begin,ii_end
93        CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown)
94        CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown)
95      ENDDO   
96
97      DO j=jj_begin,jj_end
98        CALL request_add_point(ind,ii_begin-1,j,req_dissip,right)
99        CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown)
100      ENDDO   
101
102      DO i=ii_begin+1,ii_end-1
103        CALL request_add_point(ind,i,jj_begin,req_dissip,right)
104        CALL request_add_point(ind,i,jj_end,req_dissip,right)
105      ENDDO
106   
107      DO j=jj_begin+1,jj_end-1
108        CALL request_add_point(ind,ii_begin,j,req_dissip,rup)
109        CALL request_add_point(ind,ii_end,j,req_dissip,rup)
110      ENDDO   
111
112      CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left)
113      CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown)
114      CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left)
115      CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown)
116     
117    ENDDO
118
119
120    cgraddiv=-1
121    cdivgrad=-1
122    cgradrot=-1
123 
124    CALL RANDOM_SEED()
125   
126    DO ind=1,ndomain
127      CALL swap_dimensions(ind)
128      CALL swap_geometry(ind)
129      u=f_u(ind)
130
131      DO j=jj_begin,jj_end
132        DO i=ii_begin,ii_end
133          n=(j-1)*iim+i
134          CALL RANDOM_NUMBER(r)
135          u(n+u_right)=r-0.5
136          CALL RANDOM_NUMBER(r)
137          u(n+u_lup)=r-0.5
138          CALL RANDOM_NUMBER(r)
139          u(n+u_ldown)=r-0.5
140       ENDDO
141      ENDDO
142       
143    ENDDO
144 
145
146
147    DO it=1,20
148
149      CALL transfert_request(f_u,req_dissip)
150      CALL transfert_request(f_u,req_dissip)
151
152      dumax=0
153      DO ind=1,ndomain
154        CALL swap_dimensions(ind)
155        CALL swap_geometry(ind)
156        u=f_u(ind)
157        du=f_du(ind)
158        CALL gradiv(u,du,1)
159      ENDDO
160      CALL transfert_request(f_du,req_dissip)
161      CALL transfert_request(f_du,req_dissip)
162     
163      DO ind=1,ndomain
164        CALL swap_dimensions(ind)
165        CALL swap_geometry(ind)
166        u=f_u(ind)
167        du=f_du(ind)
168       
169        DO j=jj_begin,jj_end
170          DO i=ii_begin,ii_end
171            n=(j-1)*iim+i
172            if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
173            if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
174            if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
175          ENDDO
176        ENDDO
177      ENDDO
178
179      DO ind=1,ndomain
180        CALL swap_dimensions(ind)
181        CALL swap_geometry(ind)
182        u=f_u(ind)
183        du=f_du(ind)
184        u=du/dumax
185      ENDDO
186      PRINT *,"gradiv : it :",it ,": dumax",dumax
187
188    ENDDO 
189    PRINT *,"gradiv : dumax",dumax
190   
191    cgraddiv=dumax**(-1./nitergdiv)
192    PRINT *,"cgraddiv : ",cgraddiv
193
194    DO ind=1,ndomain
195      CALL swap_dimensions(ind)
196      CALL swap_geometry(ind)
197      u=f_u(ind)
198
199     DO j=jj_begin,jj_end
200        DO i=ii_begin,ii_end
201          n=(j-1)*iim+i
202          CALL RANDOM_NUMBER(r)
203          u(n+u_right)=r-0.5
204          CALL RANDOM_NUMBER(r)
205          u(n+u_lup)=r-0.5
206          CALL RANDOM_NUMBER(r)
207          u(n+u_ldown)=r-0.5
208       ENDDO
209      ENDDO
210       
211    ENDDO
212
213
214    DO it=1,20
215 
216      CALL transfert_request(f_u,req_dissip)
217      CALL transfert_request(f_u,req_dissip)
218
219      dumax=0
220      DO ind=1,ndomain
221        CALL swap_dimensions(ind)
222        CALL swap_geometry(ind)
223        u=f_u(ind)
224        du=f_du(ind)
225        CALL gradrot(u,du,1)
226      ENDDO
227
228      CALL transfert_request(f_du,req_dissip)
229      CALL transfert_request(f_du,req_dissip)
230     
231      DO ind=1,ndomain
232        CALL swap_dimensions(ind)
233        CALL swap_geometry(ind)
234        u=f_u(ind)
235        du=f_du(ind)
236       
237        DO j=jj_begin,jj_end
238          DO i=ii_begin,ii_end
239            n=(j-1)*iim+i
240            if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
241            if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
242            if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
243          ENDDO
244        ENDDO
245      ENDDO
246     
247      DO ind=1,ndomain
248        CALL swap_dimensions(ind)
249        CALL swap_geometry(ind)
250        u=f_u(ind)
251        du=f_du(ind)
252        u=du/dumax
253      ENDDO
254     
255!      PRINT *,"gradrot : it :",it ,": dumax",(dumax+dumaxm1)/2
256      PRINT *,"gradrot : it :",it ,": dumax",dumax
257
258    ENDDO 
259    PRINT *,"gradrot : dumax",dumax
260 
261    cgradrot=dumax**(-1/nitergrot) 
262    PRINT *,"cgradrot : ",cgradrot
263   
264
265
266    DO ind=1,ndomain
267      CALL swap_dimensions(ind)
268      CALL swap_geometry(ind)
269      theta=f_theta(ind)
270 
271       DO j=jj_begin,jj_end
272          DO i=ii_begin,ii_end
273            n=(j-1)*iim+i
274            CALL RANDOM_NUMBER(r)
275            theta(n)=r-0.5
276         ENDDO
277        ENDDO
278       
279    ENDDO
280
281    DO it=1,20
282
283      CALL transfert_request(f_theta,req_i1)
284
285      dthetamax=0
286      DO ind=1,ndomain
287        CALL swap_dimensions(ind)
288        CALL swap_geometry(ind)
289        theta=f_theta(ind)
290        dtheta=f_dtheta(ind)
291        CALL divgrad(theta,dtheta,1)
292      ENDDO
293!      CALL writefield("divgrad",f_dtheta)
294
295      CALL transfert_request(f_dtheta,req_i1)
296     
297      DO ind=1,ndomain
298        CALL swap_dimensions(ind)
299        CALL swap_geometry(ind)
300        theta=f_theta(ind)
301        dtheta=f_dtheta(ind)
302       
303        DO j=jj_begin,jj_end
304          DO i=ii_begin,ii_end
305            n=(j-1)*iim+i
306            dthetamax=MAX(dthetamax,ABS(dtheta(n)))
307          ENDDO
308        ENDDO
309      ENDDO
310
311      PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
312
313      DO ind=1,ndomain
314        CALL swap_dimensions(ind)
315        CALL swap_geometry(ind)
316        theta=f_theta(ind)
317        dtheta=f_dtheta(ind)
318        theta=dtheta/dthetamax
319      ENDDO
320    ENDDO 
321
322!    CALL writefield("divgrad",f_dtheta)
323    PRINT *,"divgrad : divgrad",dthetamax
324 
325    cdivgrad=dthetamax**(-1/niterdivgrad) 
326    PRINT *,"cdivgrad : ",cdivgrad
327   
328
329
330
331
332
333     
334    fact=2
335    DO l=1,llm
336       zz= 1. - preff/presnivs(l)
337       zvert=fact-(fact-1)/(1+zz*zz)
338       tau_graddiv(l) = tau_graddiv(l)/zvert
339       tau_gradrot(l) = tau_gradrot(l)/zvert
340       tau_divgrad(l) = tau_divgrad(l)/zvert
341    ENDDO
342       
343!    idissip=INT(MIN(tetagdiv,tetagrot)/(2.*dt))
344!    idissip=MAX(1,idissip)
345!    dtdissip=idissip*dt
346!    PRINT *,"idissip",idissip," dtdissip ",dtdissip
347
348  END SUBROUTINE init_dissip 
349 
350 
351  SUBROUTINE dissip(f_ue,f_due,f_ps,f_theta_rhodz,f_dtheta_rhodz)
352  USE icosa
353  USE theta2theta_rhodz_mod
354  IMPLICIT NONE
355    TYPE(t_field),POINTER :: f_ue(:)
356    TYPE(t_field),POINTER :: f_due(:)
357    TYPE(t_field),POINTER :: f_ps(:)
358    TYPE(t_field),POINTER :: f_theta_rhodz(:)
359    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
360    REAL(rstd),POINTER         :: ue(:,:)
361    REAL(rstd),POINTER         :: due(:,:)
362    REAL(rstd),POINTER         :: ps(:)
363    REAL(rstd),POINTER         :: theta_rhodz(:,:)
364    REAL(rstd),POINTER         :: dtheta_rhodz(:,:)
365
366    REAL(rstd)                 :: theta(iim*jjm,llm)
367    REAL(rstd)                 :: du_dissip(3*iim*jjm,llm)
368    REAL(rstd)                 :: dtheta_dissip(iim*jjm,llm)
369    REAL(rstd)                 :: dtheta_rhodz_dissip(iim*jjm,llm)
370
371    INTEGER :: ind
372    INTEGER :: l,i,j,n
373   
374    CALL transfert_request(f_ue,req_e1)
375    CALL transfert_request(f_ue,req_e1)
376    CALL transfert_request(f_theta_rhodz,req_i1)
377    CALL transfert_request(f_ps,req_i1)
378
379    DO ind=1,ndomain
380      CALL swap_dimensions(ind)
381      CALL swap_geometry(ind)
382      ue=f_ue(ind)
383      due=f_due(ind) 
384      ps=f_ps(ind) 
385      theta_rhodz=f_theta_rhodz(ind)
386      dtheta_rhodz=f_dtheta_rhodz(ind) 
387
388!$OMP PARALLEL DEFAULT(SHARED)
389     CALL compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz)
390!$OMP END PARALLEL
391     
392    ENDDO
393     
394  END SUBROUTINE dissip
395
396
397  SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz)
398  USE icosa
399  USE theta2theta_rhodz_mod
400  IMPLICIT NONE
401    REAL(rstd)         :: ue(3*iim*jjm,llm)
402    REAL(rstd)         :: due(3*iim*jjm,llm)
403    REAL(rstd)         :: ps(iim*jjm)
404    REAL(rstd)         :: theta_rhodz(iim*jjm,llm)
405    REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm)
406
407    REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:)
408    REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:)
409    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:)
410    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:)
411
412    INTEGER :: ind
413    INTEGER :: l,i,j,n
414
415!$OMP BARRIER
416!$OMP MASTER
417      ALLOCATE(theta(iim*jjm,llm))
418      ALLOCATE(du_dissip(3*iim*jjm,llm))
419      ALLOCATE(dtheta_dissip(iim*jjm,llm))
420      ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm))
421!$OMP END MASTER
422!$OMP BARRIER
423   
424      CALL gradiv(ue,du_dissip,llm)
425      DO l=1,llm
426!$OMP DO
427        DO j=jj_begin,jj_end
428          DO i=ii_begin,ii_end
429            n=(j-1)*iim+i
430            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5
431            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5
432            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5
433          ENDDO
434        ENDDO
435      ENDDO
436     
437      CALL gradrot(ue,du_dissip,llm)
438
439      DO l=1,llm
440!$OMP DO
441        DO j=jj_begin,jj_end
442          DO i=ii_begin,ii_end
443            n=(j-1)*iim+i
444            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5
445            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5
446            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5
447          ENDDO
448        ENDDO
449      ENDDO 
450     
451      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
452      CALL divgrad(theta,dtheta_dissip,llm)
453      CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0)
454
455      DO l=1,llm
456!$OMP DO
457        DO j=jj_begin,jj_end
458          DO i=ii_begin,ii_end
459            n=(j-1)*iim+i
460            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5
461          ENDDO
462        ENDDO
463      ENDDO 
464
465!$OMP BARRIER
466!$OMP MASTER
467      DEALLOCATE(theta)
468      DEALLOCATE(du_dissip)
469      DEALLOCATE(dtheta_dissip)
470      DEALLOCATE(dtheta_rhodz_dissip)
471!$OMP END MASTER
472!$OMP BARRIER     
473
474  END SUBROUTINE compute_dissip
475 
476 
477  SUBROUTINE gradiv(ue,gradivu_e,ll)
478  USE icosa
479  IMPLICIT NONE
480    INTEGER,INTENT(IN)     :: ll
481    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
482    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,ll)
483    REAL(rstd),SAVE,ALLOCATABLE :: divu_i(:,:)
484   
485    INTEGER :: i,j,n,l
486
487!$OMP BARRIER
488!$OMP MASTER
489      ALLOCATE(divu_i(iim*jjm,ll))
490!$OMP END MASTER
491!$OMP BARRIER
492
493    DO l=1,ll
494!$OMP DO
495      DO j=jj_begin,jj_end
496        DO i=ii_begin,ii_end
497          n=(j-1)*iim+i
498          divu_i(n,l)=1./Ai(n)*(ne(n,right)*ue(n+u_right,l)*le(n+u_right)  +  &
499                             ne(n,rup)*ue(n+u_rup,l)*le(n+u_rup)        +  & 
500                             ne(n,lup)*ue(n+u_lup,l)*le(n+u_lup)        +  & 
501                             ne(n,left)*ue(n+u_left,l)*le(n+u_left)     +  & 
502                             ne(n,ldown)*ue(n+u_ldown,l)*le(n+u_ldown)  +  & 
503                             ne(n,rdown)*ue(n+u_rdown,l)*le(n+u_rdown))
504        ENDDO
505      ENDDO
506    ENDDO
507   
508    DO l=1,ll
509!$OMP DO
510      DO j=jj_begin,jj_end
511        DO i=ii_begin,ii_end
512 
513          n=(j-1)*iim+i
514 
515          gradivu_e(n+u_right,l)=-1/de(n+u_right)*(ne(n,right)*divu_i(n,l)+ ne(n+t_right,left)*divu_i(n+t_right,l) )       
516
517          gradivu_e(n+u_lup,l)=-1/de(n+u_lup)*(ne(n,lup)*divu_i(n,l)+ ne(n+t_lup,rdown)*divu_i(n+t_lup,l))       
518   
519          gradivu_e(n+u_ldown,l)=-1/de(n+u_ldown)*(ne(n,ldown)*divu_i(n,l)+ne(n+t_ldown,rup)*divu_i(n+t_ldown,l) )
520
521        ENDDO
522      ENDDO
523    ENDDO
524
525    DO l=1,ll
526!$OMP DO
527      DO j=jj_begin,jj_end
528        DO i=ii_begin,ii_end
529          n=(j-1)*iim+i
530          gradivu_e(n+u_right,l)=gradivu_e(n+u_right,l)*cgraddiv
531          gradivu_e(n+u_lup,l)=gradivu_e(n+u_lup,l)*cgraddiv
532          gradivu_e(n+u_ldown,l)=gradivu_e(n+u_ldown,l)*cgraddiv
533        ENDDO
534      ENDDO
535    ENDDO
536
537!$OMP BARRIER
538!$OMP MASTER
539      DEALLOCATE(divu_i)
540!$OMP END MASTER
541!$OMP BARRIER
542
543   
544  END SUBROUTINE gradiv
545 
546  SUBROUTINE divgrad(theta,divgrad_i,ll)
547  USE icosa
548  IMPLICIT NONE
549    INTEGER,INTENT(IN)     :: ll
550    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,ll)
551    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,ll)
552    REAL(rstd),SAVE,ALLOCATABLE :: grad_e(:,:)
553
554    INTEGER :: i,j,n,l
555
556!$OMP BARRIER
557!$OMP MASTER
558      ALLOCATE(grad_e(3*iim*jjm,ll))
559!$OMP END MASTER
560!$OMP BARRIER
561       
562    DO l=1,ll
563!$OMP DO
564      DO j=jj_begin-1,jj_end+1
565        DO i=ii_begin-1,ii_end+1
566   
567          n=(j-1)*iim+i
568 
569          grad_e(n+u_right,l)=-1/de(n+u_right)*(ne(n,right)*theta(n,l)+ ne(n+t_right,left)*theta(n+t_right,l) )       
570
571          grad_e(n+u_lup,l)=-1/de(n+u_lup)*(ne(n,lup)*theta(n,l)+ ne(n+t_lup,rdown)*theta(n+t_lup,l ))       
572   
573          grad_e(n+u_ldown,l)=-1/de(n+u_ldown)*(ne(n,ldown)*theta(n,l)+ne(n+t_ldown,rup)*theta(n+t_ldown,l) )
574
575        ENDDO
576      ENDDO
577    ENDDO
578   
579   
580    DO l=1,ll
581!$OMP DO
582      DO j=jj_begin,jj_end
583        DO i=ii_begin,ii_end
584          n=(j-1)*iim+i
585          divgrad_i(n,l)=1./Ai(n)*(ne(n,right)*grad_e(n+u_right,l)*le(n+u_right)  +  &
586                             ne(n,rup)*grad_e(n+u_rup,l)*le(n+u_rup)              +  & 
587                             ne(n,lup)*grad_e(n+u_lup,l)*le(n+u_lup)              +  & 
588                             ne(n,left)*grad_e(n+u_left,l)*le(n+u_left)           +  & 
589                             ne(n,ldown)*grad_e(n+u_ldown,l)*le(n+u_ldown)        +  & 
590                             ne(n,rdown)*grad_e(n+u_rdown,l)*le(n+u_rdown))
591        ENDDO
592      ENDDO
593    ENDDO
594   
595    DO l=1,ll
596!$OMP DO
597      DO j=jj_begin,jj_end
598        DO i=ii_begin,ii_end
599          n=(j-1)*iim+i
600          divgrad_i(n,l)=divgrad_i(n,l)*cdivgrad
601        ENDDO
602      ENDDO
603    ENDDO
604
605!$OMP BARRIER
606!$OMP MASTER
607      DEALLOCATE(grad_e)
608!$OMP END MASTER
609!$OMP BARRIER
610
611  END SUBROUTINE divgrad
612
613   
614  SUBROUTINE gradrot(ue,gradrot_e,ll)
615  USE icosa
616  IMPLICIT NONE
617    INTEGER,INTENT(IN)     :: ll
618    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
619    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,ll)
620    REAL(rstd),SAVE,ALLOCATABLE  :: rot_v(:,:)
621
622    INTEGER :: i,j,n,l
623
624!$OMP BARRIER
625!$OMP MASTER
626      ALLOCATE(rot_v(2*iim*jjm,ll))
627!$OMP END MASTER
628!$OMP BARRIER
629     
630    DO l=1,ll
631!$OMP DO
632      DO j=jj_begin-1,jj_end+1
633        DO i=ii_begin-1,ii_end+1
634          n=(j-1)*iim+i
635       
636          rot_v(n+z_up,l)= 1./Av(n+z_up)*(  ne(n,rup)*ue(n+u_rup,l)*de(n+u_rup)                     &
637                                + ne(n+t_rup,left)*ue(n+t_rup+u_left,l)*de(n+t_rup+u_left)          &
638                                - ne(n,lup)*ue(n+u_lup,l)*de(n+u_lup) ) 
639                             
640          rot_v(n+z_down,l) = 1./Av(n+z_down)*( ne(n,ldown)*ue(n+u_ldown,l)*de(n+u_ldown)                 &
641                                     + ne(n+t_ldown,right)*ue(n+t_ldown+u_right,l)*de(n+t_ldown+u_right)  &
642                                     - ne(n,rdown)*ue(n+u_rdown,l)*de(n+u_rdown) )
643 
644        ENDDO
645      ENDDO
646    ENDDO                             
647 
648    DO l=1,ll
649!$OMP DO
650      DO j=jj_begin,jj_end
651        DO i=ii_begin,ii_end
652          n=(j-1)*iim+i
653 
654          gradrot_e(n+u_right,l)=1/le(n+u_right)*ne(n,right)*(rot_v(n+z_rdown,l)-rot_v(n+z_rup,l)) 
655         
656          gradrot_e(n+u_lup,l)=1/le(n+u_lup)*ne(n,lup)*(rot_v(n+z_up,l)-rot_v(n+z_lup,l)) 
657       
658          gradrot_e(n+u_ldown,l)=1/le(n+u_ldown)*ne(n,ldown)*(rot_v(n+z_ldown,l)-rot_v(n+z_down,l))
659       
660        ENDDO
661      ENDDO
662    ENDDO
663
664    DO l=1,ll
665!$OMP DO
666      DO j=jj_begin,jj_end
667        DO i=ii_begin,ii_end
668          n=(j-1)*iim+i
669   
670          gradrot_e(n+u_right,l)=gradrot_e(n+u_right,l)*cgradrot       
671          gradrot_e(n+u_lup,l)=gradrot_e(n+u_lup,l)*cgradrot
672          gradrot_e(n+u_ldown,l)=gradrot_e(n+u_ldown,l)*cgradrot
673       
674        ENDDO
675      ENDDO
676    ENDDO 
677   
678!$OMP BARRIER
679!$OMP MASTER
680      DEALLOCATE(rot_v)
681!$OMP END MASTER
682!$OMP BARRIER
683   
684  END SUBROUTINE gradrot
685
686
687END MODULE dissip_gcm_mod
688       
Note: See TracBrowser for help on using the repository browser.