source: trunk/LMDZ.COMMON/libf/dyn3d/integrd.F @ 1422

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

File size: 7.1 KB
Line 
1!
2! $Id: integrd.F 1616 2012-02-17 11:59:00Z emillour $
3!
4      SUBROUTINE integrd
5     $  (  nq,vcovm1,ucovm1,tetam1,psm1,massem1,
6     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis !,finvmaold
7     &  )
8
9      use control_mod, only : planet_type
10      USE comvert_mod, ONLY: ap,bp
11      USE comconst_mod, ONLY: pi
12      USE logic_mod, ONLY: leapf
13      USE temps_mod, ONLY: dt
14
15      IMPLICIT NONE
16
17
18c=======================================================================
19c
20c   Auteur:  P. Le Van
21c   -------
22c
23c   objet:
24c   ------
25c
26c   Incrementation des tendances dynamiques
27c
28c=======================================================================
29c-----------------------------------------------------------------------
30c   Declarations:
31c   -------------
32
33#include "dimensions.h"
34#include "paramet.h"
35#include "comgeom.h"
36#include "iniprint.h"
37
38c   Arguments:
39c   ----------
40
41      integer,intent(in) :: nq ! number of tracers to handle in this routine
42      real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
43      real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
44      real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
45      real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
46      real,intent(inout) :: ps(ip1jmp1) ! surface pressure
47      real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
48      real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
49      ! values at previous time step
50      real,intent(inout) :: vcovm1(ip1jm,llm)
51      real,intent(inout) :: ucovm1(ip1jmp1,llm)
52      real,intent(inout) :: tetam1(ip1jmp1,llm)
53      real,intent(inout) :: psm1(ip1jmp1)
54      real,intent(inout) :: massem1(ip1jmp1,llm)
55      ! the tendencies to add
56      real,intent(in) :: dv(ip1jm,llm)
57      real,intent(in) :: du(ip1jmp1,llm)
58      real,intent(in) :: dteta(ip1jmp1,llm)
59      real,intent(in) :: dp(ip1jmp1)
60      real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
61!      real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
62
63c   Local:
64c   ------
65
66      REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
67      REAL massescr( ip1jmp1,llm )
68!      REAL finvmasse(ip1jmp1,llm)
69      REAL p(ip1jmp1,llmp1)
70      REAL tpn,tps,tppn(iim),tpps(iim)
71      REAL qpn,qps,qppn(iim),qpps(iim)
72      REAL deltap( ip1jmp1,llm )
73
74      INTEGER  l,ij,iq,i,j
75
76      REAL SSUM
77
78c-----------------------------------------------------------------------
79
80      DO  l = 1,llm
81        DO  ij = 1,iip1
82         ucov(    ij    , l) = 0.
83         ucov( ij +ip1jm, l) = 0.
84         uscr(     ij      ) = 0.
85         uscr( ij +ip1jm   ) = 0.
86        ENDDO
87      ENDDO
88
89
90c    ............    integration  de       ps         ..............
91
92      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
93
94      DO ij = 1,ip1jmp1
95       pscr (ij)    = ps(ij)
96       ps (ij)      = psm1(ij) + dt * dp(ij)
97      ENDDO
98c
99      DO ij = 1,ip1jmp1
100        IF( ps(ij).LT.0. ) THEN
101         write(lunout,*) "integrd: negative surface pressure ",ps(ij)
102         write(lunout,*) " at node ij =", ij
103         ! since ij=j+(i-1)*jjp1 , we have
104         j=modulo(ij,jjp1)
105         i=1+(ij-j)/jjp1
106         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
107     &                   " lat = ",rlatu(j)*180./pi, " deg"
108         write(lunout,*) " psm1(ij)=",psm1(ij)," dt=",dt,
109     &                   " dp(ij)=",dp(ij)
110         call abort_gcm("integrd", "", 1)
111        ENDIF
112      ENDDO
113c
114      DO  ij    = 1, iim
115       tppn(ij) = aire(   ij   ) * ps(  ij    )
116       tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
117      ENDDO
118       tpn      = SSUM(iim,tppn,1)/apoln
119       tps      = SSUM(iim,tpps,1)/apols
120      DO ij   = 1, iip1
121       ps(   ij   )  = tpn
122       ps(ij+ip1jm)  = tps
123      ENDDO
124c
125c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
126c
127      CALL pression ( ip1jmp1, ap, bp, ps, p )
128      CALL massdair (     p  , masse         )
129
130! Ehouarn : we don't use/need finvmaold and finvmasse,
131!           so might as well not compute them
132!      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
133!      CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1  )
134c
135
136c    ............   integration  de  ucov, vcov,  h     ..............
137
138      DO l = 1,llm
139
140       DO ij = iip2,ip1jm
141        uscr( ij )   =  ucov( ij,l )
142        ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
143       ENDDO
144
145       DO ij = 1,ip1jm
146        vscr( ij )   =  vcov( ij,l )
147        vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
148       ENDDO
149
150       DO ij = 1,ip1jmp1
151        hscr( ij )    =  teta(ij,l)
152        teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
153     &                + dt * dteta(ij,l) / masse(ij,l)
154       ENDDO
155
156c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
157c
158c
159       DO  ij   = 1, iim
160        tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
161        tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
162       ENDDO
163        tpn      = SSUM(iim,tppn,1)/apoln
164        tps      = SSUM(iim,tpps,1)/apols
165
166       DO ij   = 1, iip1
167        teta(   ij   ,l)  = tpn
168        teta(ij+ip1jm,l)  = tps
169       ENDDO
170c
171
172       IF(leapf)  THEN
173         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
174         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
175         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
176       END IF
177
178      ENDDO ! of DO l = 1,llm
179
180
181c
182c   .......  integration de   q   ......
183c
184c$$$      IF( iadv(1).NE.3.AND.iadv(2).NE.3 )    THEN
185c$$$c
186c$$$       IF( forward. OR . leapf )  THEN
187c$$$        DO iq = 1,2
188c$$$        DO  l = 1,llm
189c$$$        DO ij = 1,ip1jmp1
190c$$$        q(ij,l,iq) = ( q(ij,l,iq)*finvmaold(ij,l) + dtvr *dq(ij,l,iq) )/
191c$$$     $                            finvmasse(ij,l)
192c$$$        ENDDO
193c$$$        ENDDO
194c$$$        ENDDO
195c$$$       ELSE
196c$$$         DO iq = 1,2
197c$$$         DO  l = 1,llm
198c$$$         DO ij = 1,ip1jmp1
199c$$$         q( ij,l,iq ) = q( ij,l,iq ) * finvmaold(ij,l) / finvmasse(ij,l)
200c$$$         ENDDO
201c$$$         ENDDO
202c$$$         ENDDO
203c$$$
204c$$$       END IF
205c$$$c
206c$$$      ENDIF
207
208      if (planet_type.eq."earth") then
209! Earth-specific treatment of first 2 tracers (water)
210        DO l = 1, llm
211          DO ij = 1, ip1jmp1
212            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
213          ENDDO
214        ENDDO
215
216        CALL qminimum( q, nq, deltap )
217
218c
219c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
220c
221
222       DO iq = 1, nq
223        DO l = 1, llm
224
225           DO ij = 1, iim
226             qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
227             qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
228           ENDDO
229             qpn  =  SSUM(iim,qppn,1)/apoln
230             qps  =  SSUM(iim,qpps,1)/apols
231
232           DO ij = 1, iip1
233             q(   ij   ,l,iq)  = qpn
234             q(ij+ip1jm,l,iq)  = qps
235           ENDDO
236
237        ENDDO
238       ENDDO
239
240! Ehouarn: forget about finvmaold
241!      CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
242
243      endif ! of if (planet_type.eq."earth")
244c
245c
246c     .....   FIN  de l'integration  de   q    .......
247
248c    .................................................................
249
250
251      IF( leapf )  THEN
252         CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
253         CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
254      END IF
255
256      RETURN
257      END
Note: See TracBrowser for help on using the repository browser.