source: LMDZ6/branches/LMDZ-INCA-Dyn/libf/dyn3dmem/advect_new_loc.F @ 5420

Last change on this file since 5420 was 2622, checked in by Ehouarn Millour, 8 years ago

Some code tidying: turn ener.h into ener_mod.F90
EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 7.3 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE advect_new_loc(ucov,vcov,teta,w,massebx,masseby,
5     &                        du,dv,dteta)
6      USE parallel_lmdz
7      USE write_field_loc
8      USE advect_new_mod
9      USE comconst_mod, ONLY: daysec
10      USE logic_mod, ONLY: conser
11     
12      IMPLICIT NONE
13c=======================================================================
14c
15c   Auteurs:  P. Le Van , Fr. Hourdin  .
16c   -------
17c
18c   Objet:
19c   ------
20c
21c   *************************************************************
22c   .... calcul des termes d'advection vertic.pour u,v,teta,q ...
23c   *************************************************************
24c        ces termes sont ajoutes a du,dv,dteta et dq .
25c  Modif F.Forget 03/94 : on retire q de advect
26c
27c=======================================================================
28c-----------------------------------------------------------------------
29c   Declarations:
30c   -------------
31
32      include "dimensions.h"
33      include "paramet.h"
34      include "comgeom.h"
35
36c   Arguments:
37c   ----------
38
39      REAL vcov(ijb_v:ije_v,llm),ucov(ijb_u:ije_u,llm)
40      REAL teta(ijb_u:ije_u,llm)
41      REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
42      REAL w(ijb_u:ije_u,llm)
43      REAL dv(ijb_v:ije_v,llm),du(ijb_u:ije_u,llm)
44      REAL dteta(ijb_u:ije_u,llm)
45c   Local:
46c   ------
47
48      REAL wsur2(ijb_u:ije_u)
49      REAL unsaire2(ijb_u:ije_u), ge(ijb_u:ije_u)
50      REAL deuxjour, ww, gt, uu, vv
51
52      INTEGER  ij,l,ijb,ije
53      EXTERNAL  SSUM
54      REAL      SSUM
55
56
57
58c-----------------------------------------------------------------------
59c   2. Calculs preliminaires:
60c   -------------------------
61     
62      IF (conser.AND.1==0)  THEN
63         deuxjour = 2. * daysec
64
65         DO   1  ij   = 1, ip1jmp1
66         unsaire2(ij) = unsaire(ij) * unsaire(ij)
67   1     CONTINUE
68      END IF
69
70
71c------------------  -yy ----------------------------------------------
72c   .  Calcul de     u
73
74c$OMP MASTER
75      ijb=ij_begin
76      ije=ij_end
77      if (pole_nord) ijb=ijb+iip1
78      if (pole_sud)  ije=ije-iip1
79
80      DO ij=ijb,ije
81        du2(ij,1)=0.
82        du1(ij,llm)=0.
83      ENDDO
84     
85      ijb=ij_begin
86      ije=ij_end
87      if (pole_sud)  ije=ij_end-iip1
88     
89      DO ij=ijb,ije
90        dv2(ij,1)=0.
91        dv1(ij,llm)=0.
92      ENDDO
93     
94      ijb=ij_begin
95      ije=ij_end
96
97      DO ij=ijb,ije
98        dteta2(ij,1)=0.
99        dteta1(ij,llm)=0.
100      ENDDO
101c$OMP END MASTER
102c$OMP BARRIER
103 
104c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
105      DO  l=1,llm
106         
107         ijb=ij_begin
108         ije=ij_end
109         if (pole_nord) ijb=ijb+iip1
110         if (pole_sud)  ije=ije-iip1
111         
112c         DO    ij     = iip2, ip1jmp1
113c            uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
114c         ENDDO
115
116c         DO    ij     = iip2, ip1jm
117c            uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
118c         ENDDO
119         
120         DO    ij     = ijb, ije
121                 
122           uav(ij,l)=0.25*(ucov(ij,l)+ucov(ij-iip1,l))
123     .               +0.25*(ucov(ij+iip1,l)+ucov(ij,l))
124         ENDDO
125         
126         if (pole_nord) then
127           DO      ij         = 1, iip1
128              uav(ij      ,l) = 0.
129           ENDDO
130         endif
131         
132         if (pole_sud) then
133           DO      ij         = 1, iip1
134              uav(ip1jm+ij,l) = 0.
135           ENDDO
136         endif
137         
138      ENDDO
139c$OMP END DO     
140c      call write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/)))
141     
142c------------------  -xx ----------------------------------------------
143c   .  Calcul de     v
144     
145      ijb=ij_begin
146      ije=ij_end
147      if (pole_sud)  ije=ij_end-iip1
148
149c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
150      DO  l=1,llm
151         
152         DO    ij   = ijb+1, ije
153           vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) )
154         ENDDO
155         
156         DO    ij   = ijb,ije,iip1
157          vav(ij,l) = vav(ij+iim,l)
158         ENDDO
159         
160         
161         DO    ij   = ijb, ije-1
162          vav(ij,l) = vav(ij,l) + vav(ij+1,l)
163         ENDDO
164         
165         DO    ij       = ijb, ije, iip1
166          vav(ij+iim,l) = vav(ij,l)
167         ENDDO
168         
169      ENDDO
170c$OMP END DO
171c       call write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))
172
173c-----------------------------------------------------------------------
174c$OMP BARRIER
175
176c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
177      DO 20 l = 1, llmm1
178
179
180c       ......   calcul de  - w/2.    au niveau  l+1   .......
181      ijb=ij_begin
182      ije=ij_end+iip1
183      if (pole_sud)  ije=ij_end
184     
185      DO 5   ij   = ijb, ije
186      wsur2( ij ) = - 0.5 * w( ij,l+1 )
187   5  CONTINUE
188
189
190c     .....................     calcul pour  du     ..................
191     
192      ijb=ij_begin
193      ije=ij_end
194      if (pole_nord) ijb=ijb+iip1
195      if (pole_sud)  ije=ije-iip1
196         
197      DO 6 ij = ijb ,ije-1
198      ww        = wsur2 (  ij  )     + wsur2( ij+1 )
199      uu        = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) )
200      du1(ij,l)  =  ww * ( uu - uav(ij, l ) )/massebx(ij, l )
201      du2(ij,l+1)=  ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1)
202   6  CONTINUE
203
204c     .................    calcul pour   dv      .....................
205      ijb=ij_begin
206      ije=ij_end
207      if (pole_sud)  ije=ij_end-iip1
208     
209      DO 8 ij = ijb, ije
210      ww        = wsur2( ij+iip1 )   + wsur2( ij )
211      vv        = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) )
212      dv1(ij,l)  =  ww * (vv - vav(ij, l ) )/masseby(ij, l )
213      dv2(ij,l+1)=  ww * (vv - vav(ij,l+1) )/masseby(ij,l+1)
214   8  CONTINUE
215
216c
217
218c     ............................................................
219c     ...............    calcul pour   dh      ...................
220c     ............................................................
221
222c                       ---z
223c       calcul de  - d( teta  * w )      qu'on ajoute a   dh
224c                   ...............
225        ijb=ij_begin
226        ije=ij_end
227       
228        DO 15 ij = ijb, ije
229         ww            = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) )
230         dteta1(ij, l ) =   ww
231         dteta2(ij,l+1) =   ww
232  15    CONTINUE
233
234c ym ---> conser a voir plus tard
235
236c      IF( conser)  THEN
237c       
238c        DO 17 ij = 1,ip1jmp1
239c        ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
240c  17    CONTINUE
241c        gt       = SSUM( ip1jmp1,ge,1 )
242c        gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
243c      END IF
244
245  20  CONTINUE
246c$OMP END DO
247
248      ijb=ij_begin
249      ije=ij_end
250      if (pole_nord) ijb=ijb+iip1
251      if (pole_sud)  ije=ije-iip1
252#ifdef DEBUG_IO   
253       CALL WriteField_u('du_bis',du)     
254#endif
255c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
256      DO l=1,llm
257        DO ij=ijb,ije-1
258          du(ij,l)=du(ij,l)+du2(ij,l)-du1(ij,l)
259        ENDDO
260
261        DO   ij   = ijb+iip1-1, ije, iip1
262         du( ij, l  ) = du( ij -iim, l  )
263        ENDDO
264      ENDDO
265c$OMP END DO NOWAIT
266#ifdef DEBUG_IO   
267      CALL WriteField_u('du1',du1)     
268      CALL WriteField_u('du2',du2)     
269      CALL WriteField_u('du_bis',du)     
270#endif
271      ijb=ij_begin
272      ije=ij_end
273      if (pole_sud)  ije=ij_end-iip1
274
275c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
276      DO l=1,llm
277        DO ij=ijb,ije
278          dv(ij,l)=dv(ij,l)+dv2(ij,l)-dv1(ij,l)
279        ENDDO
280      ENDDO
281c$OMP END DO NOWAIT     
282      ijb=ij_begin
283      ije=ij_end
284
285c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
286      DO l=1,llm
287        DO ij=ijb,ije
288          dteta(ij,l)=dteta(ij,l)+dteta2(ij,l)-dteta1(ij,l)
289        ENDDO
290      ENDDO
291c$OMP END DO NOWAIT     
292
293      RETURN
294      END
Note: See TracBrowser for help on using the repository browser.