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

Last change on this file since 706 was 630, checked in by Laurent Fairhead, 20 years ago

Import d'une version parallele de la dynamique YM
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.6 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      EXTERNAL  gradiv ,nXgrarot,divgrad,initial0
52      EXTERNAL  gradiv2,nXgraro2,divgrad2,SSUM
53      integer :: ijb,ije
54c-----------------------------------------------------------------------
55c   initialisations:
56c   ----------------
57
58      DO l=1,llm
59         te1dt(l) = tetaudiv(l) * dtdiss
60         te2dt(l) = tetaurot(l) * dtdiss
61         te3dt(l) = tetah(l)    * dtdiss
62      ENDDO
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     
70      du(ijb:ije,:)=0
71      dh(ijb:ije,:)=0
72     
73      if (pole_sud) ije=ij_end-iip1
74     
75      dv(ijb:ije,:)=0
76     
77c-----------------------------------------------------------------------
78c   Calcul de la dissipation:
79c   -------------------------
80
81c   Calcul de la partie   grad  ( div ) :
82c   -------------------------------------
83     
84     
85     
86      IF(lstardis) THEN
87c      IF (.FALSE.) THEN
88         CALL gradiv2_p( llm,ucov,vcov,nitergdiv,gdx,gdy )
89      ELSE
90         CALL gradiv_p ( llm,ucov,vcov,nitergdiv,gdx,gdy )
91      ENDIF
92
93c      call write_field3d_p('gdx',reshape(gdx,(/iip1,jjp1,llm/)))
94c     call write_field3d_p('gdy',reshape(gdy,(/iip1,jjm,llm/)))
95c       stop
96
97      ijb=ij_begin
98      ije=ij_end
99      if (pole_sud) ije=ij_end-iip1
100
101      DO l=1,llm
102         if (pole_nord) then
103           DO ij = 1, iip1
104              gdx(     ij ,l) = 0.
105           ENDDO
106         endif
107         
108         if (pole_sud) then
109           DO ij = 1, iip1
110              gdx(ij+ip1jm,l) = 0.
111           ENDDO
112         endif
113         
114         if (pole_nord) ijb=ij_begin+iip1
115         DO ij = ijb,ije
116            du(ij,l) = du(ij,l) - te1dt(l) *gdx(ij,l)
117         ENDDO
118
119         if (pole_nord) ijb=ij_begin
120         DO ij = ijb,ije
121            dv(ij,l) = dv(ij,l) - te1dt(l) *gdy(ij,l)
122         ENDDO
123
124       ENDDO
125
126c   calcul de la partie   n X grad ( rot ):
127c   ---------------------------------------
128
129      IF(lstardis) THEN
130c      IF (.FALSE.) THEN
131         CALL nxgraro2_p( llm,ucov, vcov, nitergrot,grx,gry )
132      ELSE
133         CALL nxgrarot_p( llm,ucov, vcov, nitergrot,grx,gry )
134      ENDIF
135
136c      call write_field3d_p('grx',reshape(grx,(/iip1,jjp1,llm/)))
137c      call write_field3d_p('gry',reshape(gry,(/iip1,jjm,llm/)))
138c          stop
139
140
141      ijb=ij_begin
142      ije=ij_end
143      if (pole_sud) ije=ij_end-iip1
144     
145      DO l=1,llm
146         
147         if (pole_nord) then
148           DO ij = 1, iip1
149              grx(ij,l) = 0.
150           ENDDO
151         endif
152         
153         if (pole_nord) ijb=ij_begin+iip1
154         DO ij = ijb,ije
155            du(ij,l) = du(ij,l) - te2dt(l) * grx(ij,l)
156         ENDDO
157         
158         if (pole_nord) ijb=ij_begin
159         DO ij =  ijb, ije
160            dv(ij,l) = dv(ij,l) - te2dt(l) * gry(ij,l)
161         ENDDO
162     
163      ENDDO
164
165c   calcul de la partie   div ( grad ):
166c   -----------------------------------
167
168       
169      IF(lstardis) THEN
170c      IF (.FALSE.) THEN
171   
172      ijb=ij_begin
173      ije=ij_end
174     
175       DO l = 1, llm
176          DO ij = ijb, ije
177            deltapres(ij,l) = AMAX1( 0.,  p(ij,l) - p(ij,l+1) )
178          ENDDO
179       ENDDO
180
181         CALL divgrad2_p( llm,teta, deltapres  ,niterh, gdx )
182      ELSE
183         CALL divgrad_p ( llm,teta, niterh, gdx        )
184      ENDIF
185
186c      call write_field3d_p('gdx2',reshape(gdx,(/iip1,jmp1,llm/)))
187c      stop
188
189      ijb=ij_begin
190      ije=ij_end
191     
192      DO l = 1,llm
193         DO ij = ijb,ije
194            dh( ij,l ) = dh( ij,l ) - te3dt(l) * gdx( ij,l )
195         ENDDO
196      ENDDO
197
198      RETURN
199      END
Note: See TracBrowser for help on using the repository browser.