source: LMDZ4/trunk/libf/dyn3dpar/dissip_p.F @ 3635

Last change on this file since 3635 was 985, checked in by Laurent Fairhead, 16 years ago

Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.7 KB
Line 
1      SUBROUTINE dissip_p( vcov,ucov,teta,p, dv,du,dh )
2c
3      USE parallel
4      USE write_field_p
5      IMPLICIT NONE
6
7
8c ..  Avec nouveaux operateurs star :  gradiv2 , divgrad2, nxgraro2  ...
9c                                 (  10/01/98  )
10
11c=======================================================================
12c
13c   Auteur:  P. Le Van
14c   -------
15c
16c   Objet:
17c   ------
18c
19c   Dissipation horizontale
20c
21c=======================================================================
22c-----------------------------------------------------------------------
23c   Declarations:
24c   -------------
25
26#include "dimensions.h"
27#include "paramet.h"
28#include "comconst.h"
29#include "comgeom.h"
30#include "comdissnew.h"
31#include "comdissipn.h"
32
33c   Arguments:
34c   ----------
35
36      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
37      REAL  p( ip1jmp1,llmp1 )
38      REAL dv(ip1jm,llm),du(ip1jmp1,llm),dh(ip1jmp1,llm)
39
40c   Local:
41c   ------
42
43      REAL gdx(ip1jmp1,llm),gdy(ip1jm,llm)
44      REAL grx(ip1jmp1,llm),gry(ip1jm,llm)
45      REAL te1dt(llm),te2dt(llm),te3dt(llm)
46      REAL deltapres(ip1jmp1,llm)
47
48      INTEGER l,ij
49
50      REAL  SSUM
51      integer :: ijb,ije
52c-----------------------------------------------------------------------
53c   initialisations:
54c   ----------------
55
56c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
57      DO l=1,llm
58         te1dt(l) = tetaudiv(l) * dtdiss
59         te2dt(l) = tetaurot(l) * dtdiss
60         te3dt(l) = tetah(l)    * dtdiss
61      ENDDO
62c$OMP END DO NOWAIT
63c      CALL initial0( ijp1llm, du )
64c      CALL initial0( ijmllm , dv )
65c      CALL initial0( ijp1llm, dh )
66     
67      ijb=ij_begin
68      ije=ij_end
69
70c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
71      DO l=1,llm
72        du(ijb:ije,l)=0
73        dh(ijb:ije,l)=0
74      ENDDO
75c$OMP END DO NOWAIT
76     
77      if (pole_sud) ije=ij_end-iip1
78
79c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
80      DO l=1,llm
81        dv(ijb:ije,l)=0
82      ENDDO
83c$OMP END DO NOWAIT
84     
85c-----------------------------------------------------------------------
86c   Calcul de la dissipation:
87c   -------------------------
88
89c   Calcul de la partie   grad  ( div ) :
90c   -------------------------------------
91     
92     
93     
94      IF(lstardis) THEN
95c      IF (.FALSE.) THEN
96         CALL gradiv2_p( llm,ucov,vcov,nitergdiv,gdx,gdy )
97      ELSE
98         CALL gradiv_p ( llm,ucov,vcov,nitergdiv,gdx,gdy )
99      ENDIF
100
101
102      ijb=ij_begin
103      ije=ij_end
104      if (pole_sud) ije=ij_end-iip1
105
106c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
107      DO l=1,llm
108         if (pole_nord) then
109           DO ij = 1, iip1
110              gdx(     ij ,l) = 0.
111           ENDDO
112         endif
113         
114         if (pole_sud) then
115           DO ij = 1, iip1
116              gdx(ij+ip1jm,l) = 0.
117           ENDDO
118         endif
119         
120         if (pole_nord) ijb=ij_begin+iip1
121         DO ij = ijb,ije
122            du(ij,l) = du(ij,l) - te1dt(l) *gdx(ij,l)
123         ENDDO
124
125         if (pole_nord) ijb=ij_begin
126         DO ij = ijb,ije
127            dv(ij,l) = dv(ij,l) - te1dt(l) *gdy(ij,l)
128         ENDDO
129
130       ENDDO
131c$OMP END DO NOWAIT
132c   calcul de la partie   n X grad ( rot ):
133c   ---------------------------------------
134
135      IF(lstardis) THEN
136c      IF (.FALSE.) THEN
137         CALL nxgraro2_p( llm,ucov, vcov, nitergrot,grx,gry )
138      ELSE
139         CALL nxgrarot_p( llm,ucov, vcov, nitergrot,grx,gry )
140      ENDIF
141
142
143
144      ijb=ij_begin
145      ije=ij_end
146      if (pole_sud) ije=ij_end-iip1
147
148c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
149      DO l=1,llm
150         
151         if (pole_nord) then
152           DO ij = 1, iip1
153              grx(ij,l) = 0.
154           ENDDO
155         endif
156         
157         if (pole_nord) ijb=ij_begin+iip1
158         DO ij = ijb,ije
159            du(ij,l) = du(ij,l) - te2dt(l) * grx(ij,l)
160         ENDDO
161         
162         if (pole_nord) ijb=ij_begin
163         DO ij =  ijb, ije
164            dv(ij,l) = dv(ij,l) - te2dt(l) * gry(ij,l)
165         ENDDO
166     
167      ENDDO
168c$OMP END DO NOWAIT
169
170c   calcul de la partie   div ( grad ):
171c   -----------------------------------
172
173       
174      IF(lstardis) THEN
175c      IF (.FALSE.) THEN
176   
177      ijb=ij_begin
178      ije=ij_end
179
180c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
181       DO l = 1, llm
182          DO ij = ijb, ije
183            deltapres(ij,l) = AMAX1( 0.,  p(ij,l) - p(ij,l+1) )
184          ENDDO
185       ENDDO
186c$OMP END DO NOWAIT
187         CALL divgrad2_p( llm,teta, deltapres  ,niterh, gdx )
188      ELSE
189         CALL divgrad_p ( llm,teta, niterh, gdx        )
190      ENDIF
191
192c      call write_field3d_p('gdx2',reshape(gdx,(/iip1,jmp1,llm/)))
193c      stop
194
195      ijb=ij_begin
196      ije=ij_end
197     
198c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
199      DO l = 1,llm
200         DO ij = ijb,ije
201            dh( ij,l ) = dh( ij,l ) - te3dt(l) * gdx( ij,l )
202         ENDDO
203      ENDDO
204c$OMP END DO NOWAIT
205
206      RETURN
207      END
Note: See TracBrowser for help on using the repository browser.