source: LMDZ4/branches/V3_test/libf/dyn3dpar/dissip_p.F @ 3733

Last change on this file since 3733 was 709, checked in by Laurent Fairhead, 18 years ago

Nouvelles versions de la dynamique YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.0 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
101c      call write_field3d_p('gdx',reshape(gdx,(/iip1,jjp1,llm/)))
102c     call write_field3d_p('gdy',reshape(gdy,(/iip1,jjm,llm/)))
103c       stop
104
105      ijb=ij_begin
106      ije=ij_end
107      if (pole_sud) ije=ij_end-iip1
108
109c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
110      DO l=1,llm
111         if (pole_nord) then
112           DO ij = 1, iip1
113              gdx(     ij ,l) = 0.
114           ENDDO
115         endif
116         
117         if (pole_sud) then
118           DO ij = 1, iip1
119              gdx(ij+ip1jm,l) = 0.
120           ENDDO
121         endif
122         
123         if (pole_nord) ijb=ij_begin+iip1
124         DO ij = ijb,ije
125            du(ij,l) = du(ij,l) - te1dt(l) *gdx(ij,l)
126         ENDDO
127
128         if (pole_nord) ijb=ij_begin
129         DO ij = ijb,ije
130            dv(ij,l) = dv(ij,l) - te1dt(l) *gdy(ij,l)
131         ENDDO
132
133       ENDDO
134c$OMP END DO NOWAIT
135c   calcul de la partie   n X grad ( rot ):
136c   ---------------------------------------
137
138      IF(lstardis) THEN
139c      IF (.FALSE.) THEN
140         CALL nxgraro2_p( llm,ucov, vcov, nitergrot,grx,gry )
141      ELSE
142         CALL nxgrarot_p( llm,ucov, vcov, nitergrot,grx,gry )
143      ENDIF
144
145c      call write_field3d_p('grx',reshape(grx,(/iip1,jjp1,llm/)))
146c      call write_field3d_p('gry',reshape(gry,(/iip1,jjm,llm/)))
147c          stop
148
149
150      ijb=ij_begin
151      ije=ij_end
152      if (pole_sud) ije=ij_end-iip1
153
154c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
155      DO l=1,llm
156         
157         if (pole_nord) then
158           DO ij = 1, iip1
159              grx(ij,l) = 0.
160           ENDDO
161         endif
162         
163         if (pole_nord) ijb=ij_begin+iip1
164         DO ij = ijb,ije
165            du(ij,l) = du(ij,l) - te2dt(l) * grx(ij,l)
166         ENDDO
167         
168         if (pole_nord) ijb=ij_begin
169         DO ij =  ijb, ije
170            dv(ij,l) = dv(ij,l) - te2dt(l) * gry(ij,l)
171         ENDDO
172     
173      ENDDO
174c$OMP END DO NOWAIT
175
176c   calcul de la partie   div ( grad ):
177c   -----------------------------------
178
179       
180      IF(lstardis) THEN
181c      IF (.FALSE.) THEN
182   
183      ijb=ij_begin
184      ije=ij_end
185
186c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
187       DO l = 1, llm
188          DO ij = ijb, ije
189            deltapres(ij,l) = AMAX1( 0.,  p(ij,l) - p(ij,l+1) )
190          ENDDO
191       ENDDO
192c$OMP END DO NOWAIT
193         CALL divgrad2_p( llm,teta, deltapres  ,niterh, gdx )
194      ELSE
195         CALL divgrad_p ( llm,teta, niterh, gdx        )
196      ENDIF
197
198c      call write_field3d_p('gdx2',reshape(gdx,(/iip1,jmp1,llm/)))
199c      stop
200
201      ijb=ij_begin
202      ije=ij_end
203     
204c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
205      DO l = 1,llm
206         DO ij = ijb,ije
207            dh( ij,l ) = dh( ij,l ) - te3dt(l) * gdx( ij,l )
208         ENDDO
209      ENDDO
210c$OMP END DO NOWAIT
211
212      RETURN
213      END
Note: See TracBrowser for help on using the repository browser.