source: LMDZ5/trunk/libf/dyn3dpar/advect_p.F @ 3771

Last change on this file since 3771 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
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.8 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE advect_p(ucov,vcov,teta,w,massebx,masseby,du,dv,dteta)
5      USE parallel_lmdz
6      USE write_field_p
7      USE comconst_mod, ONLY: daysec
8      USE logic_mod, ONLY: conser
9     
10      IMPLICIT NONE
11c=======================================================================
12c
13c   Auteurs:  P. Le Van , Fr. Hourdin  .
14c   -------
15c
16c   Objet:
17c   ------
18c
19c   *************************************************************
20c   .... calcul des termes d'advection vertic.pour u,v,teta,q ...
21c   *************************************************************
22c        ces termes sont ajoutes a du,dv,dteta et dq .
23c  Modif F.Forget 03/94 : on retire q de advect
24c
25c=======================================================================
26c-----------------------------------------------------------------------
27c   Declarations:
28c   -------------
29
30#include "dimensions.h"
31#include "paramet.h"
32#include "comgeom.h"
33
34c   Arguments:
35c   ----------
36
37      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
38      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm),w(ip1jmp1,llm)
39      REAL dv(ip1jm,llm),du(ip1jmp1,llm),dteta(ip1jmp1,llm)
40
41c   Local:
42c   ------
43
44      REAL uav(ip1jmp1,llm),vav(ip1jm,llm),wsur2(ip1jmp1)
45      REAL unsaire2(ip1jmp1), ge(ip1jmp1)
46      REAL deuxjour, ww, gt, uu, vv
47
48      INTEGER  ij,l,ijb,ije
49
50      EXTERNAL  SSUM
51      REAL      SSUM
52
53c-----------------------------------------------------------------------
54c   2. Calculs preliminaires:
55c   -------------------------
56
57      IF (conser)  THEN
58         deuxjour = 2. * daysec
59
60         DO   1  ij   = 1, ip1jmp1
61         unsaire2(ij) = unsaire(ij) * unsaire(ij)
62   1     CONTINUE
63      END IF
64
65
66c------------------  -yy ----------------------------------------------
67c   .  Calcul de     u
68
69      DO  l=1,llm
70         
71         ijb=ij_begin
72         ije=ij_end
73         if (pole_nord) ijb=ijb+iip1
74         if (pole_sud)  ije=ije-iip1
75         
76c         DO    ij     = iip2, ip1jmp1
77c            uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
78c         ENDDO
79
80c         DO    ij     = iip2, ip1jm
81c            uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
82c         ENDDO
83         
84         DO    ij     = ijb, ije
85                 
86           uav(ij,l)=0.25*(ucov(ij,l)+ucov(ij-iip1,l))
87     .                     +0.25*(ucov(ij+iip1,l)+ucov(ij,l))
88         ENDDO
89         
90         if (pole_nord) then
91           DO      ij         = 1, iip1
92              uav(ij      ,l) = 0.
93           ENDDO
94         endif
95         
96         if (pole_sud) then
97           DO      ij         = 1, iip1
98              uav(ip1jm+ij,l) = 0.
99           ENDDO
100         endif
101         
102      ENDDO
103     
104c      call write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/)))
105     
106c------------------  -xx ----------------------------------------------
107c   .  Calcul de     v
108     
109      ijb=ij_begin
110      ije=ij_end
111      if (pole_sud)  ije=ij_end-iip1
112     
113      DO  l=1,llm
114         
115         DO    ij   = ijb+1, ije
116           vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) )
117         ENDDO
118         
119         DO    ij   = ijb,ije,iip1
120          vav(ij,l) = vav(ij+iim,l)
121         ENDDO
122         
123         
124         DO    ij   = ijb, ije-1
125          vav(ij,l) = vav(ij,l) + vav(ij+1,l)
126         ENDDO
127         
128         DO    ij       = ijb, ije, iip1
129          vav(ij+iim,l) = vav(ij,l)
130         ENDDO
131         
132      ENDDO
133c       call write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))
134c-----------------------------------------------------------------------
135
136
137     
138      DO 20 l = 1, llmm1
139
140
141c       ......   calcul de  - w/2.    au niveau  l+1   .......
142      ijb=ij_begin
143      ije=ij_end+iip1
144      if (pole_sud)  ije=ij_end
145     
146      DO 5   ij   = ijb, ije
147      wsur2( ij ) = - 0.5 * w( ij,l+1 )
148   5  CONTINUE
149
150
151c     .....................     calcul pour  du     ..................
152     
153      ijb=ij_begin
154      ije=ij_end
155      if (pole_nord) ijb=ijb+iip1
156      if (pole_sud)  ije=ije-iip1
157         
158      DO 6 ij = ijb ,ije-1
159      ww        = wsur2 (  ij  )     + wsur2( ij+1 )
160      uu        = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) )
161      du(ij,l)  = du(ij,l)   - ww * ( uu - uav(ij, l ) )/massebx(ij, l )
162      du(ij,l+1)= du(ij,l+1) + ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1)
163   6  CONTINUE
164
165c     .....  correction pour  du(iip1,j,l)  ........
166c     .....     du(iip1,j,l)= du(1,j,l)   .....
167
168CDIR$ IVDEP
169      DO   7  ij   = ijb+iip1-1, ije, iip1
170      du( ij, l  ) = du( ij -iim, l  )
171      du( ij,l+1 ) = du( ij -iim,l+1 )
172   7  CONTINUE
173
174c     .................    calcul pour   dv      .....................
175      ijb=ij_begin
176      ije=ij_end
177      if (pole_sud)  ije=ij_end-iip1
178     
179      DO 8 ij = ijb, ije
180      ww        = wsur2( ij+iip1 )   + wsur2( ij )
181      vv        = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) )
182      dv(ij,l)  = dv(ij, l ) - ww * (vv - vav(ij, l ) )/masseby(ij, l )
183      dv(ij,l+1)= dv(ij,l+1) + ww * (vv - vav(ij,l+1) )/masseby(ij,l+1)
184   8  CONTINUE
185
186c
187
188c     ............................................................
189c     ...............    calcul pour   dh      ...................
190c     ............................................................
191
192c                       ---z
193c       calcul de  - d( teta  * w )      qu'on ajoute a   dh
194c                   ...............
195        ijb=ij_begin
196        ije=ij_end
197       
198        DO 15 ij = ijb, ije
199         ww            = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) )
200         dteta(ij, l ) = dteta(ij, l )  -  ww
201         dteta(ij,l+1) = dteta(ij,l+1)  +  ww
202  15    CONTINUE
203
204c ym ---> conser a voir plus tard
205
206c      IF( conser)  THEN
207c       
208c        DO 17 ij = 1,ip1jmp1
209c        ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
210c  17    CONTINUE
211c        gt       = SSUM( ip1jmp1,ge,1 )
212c        gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
213c      END IF
214
215  20  CONTINUE
216 
217      RETURN
218      END
Note: See TracBrowser for help on using the repository browser.