source: LMDZ5/trunk/libf/dyn3d/integrd.F @ 1676

Last change on this file since 1676 was 1616, checked in by Ehouarn Millour, 13 years ago

Some cleanup around what is done during the integration step of dynamical tendencies and namely removed computation of (unused) finvmaold, thereby saving us the expense of a call to the (costly) filter at every dynamical time step.
Checked (on Vargas, in seq, omp, mpi and mixed mode) that this doesn't change the GCM results, as expected.
EM

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