source: dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3d/integrd.F

Last change on this file was 3990, checked in by millour, 8 years ago

An updated LMDZ5 (based on LMDZ rev 2786) to play with.
EM

File size: 7.1 KB
Line 
1!
2! $Id: integrd.F 2603 2016-07-25 09:31:56Z 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 comconst_mod, only: pi
11      USE logic_mod, ONLY: leapf
12      use comvert_mod, only: ap, bp
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         call abort_gcm("integrd", "", 1)
109        ENDIF
110      ENDDO
111c
112      DO  ij    = 1, iim
113       tppn(ij) = aire(   ij   ) * ps(  ij    )
114       tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
115      ENDDO
116       tpn      = SSUM(iim,tppn,1)/apoln
117       tps      = SSUM(iim,tpps,1)/apols
118      DO ij   = 1, iip1
119       ps(   ij   )  = tpn
120       ps(ij+ip1jm)  = tps
121      ENDDO
122c
123c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
124c
125      CALL pression ( ip1jmp1, ap, bp, ps, p )
126      CALL massdair (     p  , masse         )
127
128! Ehouarn : we don't use/need finvmaold and finvmasse,
129!           so might as well not compute them
130!      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
131!      CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1  )
132c
133
134c    ............   integration  de  ucov, vcov,  h     ..............
135
136      DO l = 1,llm
137
138       DO ij = iip2,ip1jm
139        uscr( ij )   =  ucov( ij,l )
140        ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
141       ENDDO
142
143       DO ij = 1,ip1jm
144        vscr( ij )   =  vcov( ij,l )
145        vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
146       ENDDO
147
148       DO ij = 1,ip1jmp1
149        hscr( ij )    =  teta(ij,l)
150        teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
151     &                + dt * dteta(ij,l) / masse(ij,l)
152       ENDDO
153
154c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
155c
156c
157       DO  ij   = 1, iim
158        tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
159        tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
160       ENDDO
161        tpn      = SSUM(iim,tppn,1)/apoln
162        tps      = SSUM(iim,tpps,1)/apols
163
164       DO ij   = 1, iip1
165        teta(   ij   ,l)  = tpn
166        teta(ij+ip1jm,l)  = tps
167       ENDDO
168c
169
170       IF(leapf)  THEN
171         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
172         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
173         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
174       END IF
175
176      ENDDO ! of DO l = 1,llm
177
178
179c
180c   .......  integration de   q   ......
181c
182c$$$      IF( iadv(1).NE.3.AND.iadv(2).NE.3 )    THEN
183c$$$c
184c$$$       IF( forward. OR . leapf )  THEN
185c$$$        DO iq = 1,2
186c$$$        DO  l = 1,llm
187c$$$        DO ij = 1,ip1jmp1
188c$$$        q(ij,l,iq) = ( q(ij,l,iq)*finvmaold(ij,l) + dtvr *dq(ij,l,iq) )/
189c$$$     $                            finvmasse(ij,l)
190c$$$        ENDDO
191c$$$        ENDDO
192c$$$        ENDDO
193c$$$       ELSE
194c$$$         DO iq = 1,2
195c$$$         DO  l = 1,llm
196c$$$         DO ij = 1,ip1jmp1
197c$$$         q( ij,l,iq ) = q( ij,l,iq ) * finvmaold(ij,l) / finvmasse(ij,l)
198c$$$         ENDDO
199c$$$         ENDDO
200c$$$         ENDDO
201c$$$
202c$$$       END IF
203c$$$c
204c$$$      ENDIF
205
206      if (planet_type.eq."earth") then
207! Earth-specific treatment of first 2 tracers (water)
208        DO l = 1, llm
209          DO ij = 1, ip1jmp1
210            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
211          ENDDO
212        ENDDO
213
214        CALL qminimum( q, nq, deltap )
215
216c
217c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
218c
219
220       DO iq = 1, nq
221        DO l = 1, llm
222
223           DO ij = 1, iim
224             qppn(ij) = aire(   ij   ) * q(   ij   ,l,iq)
225             qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
226           ENDDO
227             qpn  =  SSUM(iim,qppn,1)/apoln
228             qps  =  SSUM(iim,qpps,1)/apols
229
230           DO ij = 1, iip1
231             q(   ij   ,l,iq)  = qpn
232             q(ij+ip1jm,l,iq)  = qps
233           ENDDO
234
235        ENDDO
236       ENDDO
237
238! Ehouarn: forget about finvmaold
239!      CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
240
241      endif ! of if (planet_type.eq."earth")
242c
243c
244c     .....   FIN  de l'integration  de   q    .......
245
246c    .................................................................
247
248
249      IF( leapf )  THEN
250         CALL SCOPY (    ip1jmp1 ,  pscr   , 1,   psm1  , 1 )
251         CALL SCOPY ( ip1jmp1*llm, massescr, 1,  massem1, 1 )
252      END IF
253
254      RETURN
255      END
Note: See TracBrowser for help on using the repository browser.