source: LMDZ4/trunk/libf/dyn3d/advect.F @ 4077

Last change on this file since 4077 was 524, checked in by lmdzadmin, 20 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.5 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE advect(ucov,vcov,teta,w,massebx,masseby,du,dv,dteta)
5
6      IMPLICIT NONE
7c=======================================================================
8c
9c   Auteurs:  P. Le Van , Fr. Hourdin  .
10c   -------
11c
12c   Objet:
13c   ------
14c
15c   *************************************************************
16c   .... calcul des termes d'advection vertic.pour u,v,teta,q ...
17c   *************************************************************
18c        ces termes sont ajoutes a du,dv,dteta et dq .
19c  Modif F.Forget 03/94 : on retire q de advect
20c
21c=======================================================================
22c-----------------------------------------------------------------------
23c   Declarations:
24c   -------------
25
26#include "dimensions.h"
27#include "paramet.h"
28#include "comconst.h"
29#include "comvert.h"
30#include "comgeom.h"
31#include "logic.h"
32#include "ener.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
49
50      REAL      SSUM
51
52c-----------------------------------------------------------------------
53c   2. Calculs preliminaires:
54c   -------------------------
55
56      IF (conser)  THEN
57         deuxjour = 2. * daysec
58
59         DO   1  ij   = 1, ip1jmp1
60         unsaire2(ij) = unsaire(ij) * unsaire(ij)
61   1     CONTINUE
62      END IF
63
64
65c------------------  -yy ----------------------------------------------
66c   .  Calcul de     u
67
68      DO  l=1,llm
69         DO    ij     = iip2, ip1jmp1
70            uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
71         ENDDO
72         DO    ij     = iip2, ip1jm
73            uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
74         ENDDO
75         DO      ij         = 1, iip1
76            uav(ij      ,l) = 0.
77            uav(ip1jm+ij,l) = 0.
78         ENDDO
79      ENDDO
80
81c------------------  -xx ----------------------------------------------
82c   .  Calcul de     v
83
84      DO  l=1,llm
85         DO    ij   = 2, ip1jm
86          vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) )
87         ENDDO
88         DO    ij   = 1,ip1jm,iip1
89          vav(ij,l) = vav(ij+iim,l)
90         ENDDO
91         DO    ij   = 1, ip1jm-1
92          vav(ij,l) = vav(ij,l) + vav(ij+1,l)
93         ENDDO
94         DO    ij       = 1, ip1jm, iip1
95          vav(ij+iim,l) = vav(ij,l)
96         ENDDO
97      ENDDO
98
99c-----------------------------------------------------------------------
100
101c
102      DO 20 l = 1, llmm1
103
104
105c       ......   calcul de  - w/2.    au niveau  l+1   .......
106
107      DO 5   ij   = 1, ip1jmp1
108      wsur2( ij ) = - 0.5 * w( ij,l+1 )
109   5  CONTINUE
110
111
112c     .....................     calcul pour  du     ..................
113
114      DO 6 ij = iip2 ,ip1jm-1
115      ww        = wsur2 (  ij  )     + wsur2( ij+1 )
116      uu        = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) )
117      du(ij,l)  = du(ij,l)   - ww * ( uu - uav(ij, l ) )/massebx(ij, l )
118      du(ij,l+1)= du(ij,l+1) + ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1)
119   6  CONTINUE
120
121c     .....  correction pour  du(iip1,j,l)  ........
122c     .....     du(iip1,j,l)= du(1,j,l)   .....
123
124CDIR$ IVDEP
125      DO   7  ij   = iip1 +iip1, ip1jm, iip1
126      du( ij, l  ) = du( ij -iim, l  )
127      du( ij,l+1 ) = du( ij -iim,l+1 )
128   7  CONTINUE
129
130c     .................    calcul pour   dv      .....................
131
132      DO 8 ij = 1, ip1jm
133      ww        = wsur2( ij+iip1 )   + wsur2( ij )
134      vv        = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) )
135      dv(ij,l)  = dv(ij, l ) - ww * (vv - vav(ij, l ) )/masseby(ij, l )
136      dv(ij,l+1)= dv(ij,l+1) + ww * (vv - vav(ij,l+1) )/masseby(ij,l+1)
137   8  CONTINUE
138
139c
140
141c     ............................................................
142c     ...............    calcul pour   dh      ...................
143c     ............................................................
144
145c                       ---z
146c       calcul de  - d( teta  * w )      qu'on ajoute a   dh
147c                   ...............
148
149        DO 15 ij = 1, ip1jmp1
150         ww            = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) )
151         dteta(ij, l ) = dteta(ij, l )  -  ww
152         dteta(ij,l+1) = dteta(ij,l+1)  +  ww
153  15    CONTINUE
154
155      IF( conser)  THEN
156        DO 17 ij = 1,ip1jmp1
157        ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
158  17    CONTINUE
159        gt       = SSUM( ip1jmp1,ge,1 )
160        gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
161      END IF
162
163  20  CONTINUE
164 
165      RETURN
166      END
Note: See TracBrowser for help on using the repository browser.