/[lmdze]/trunk/Sources/dyn3d/integrd.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/integrd.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 28 by guez, Fri Mar 26 18:33:04 2010 UTC revision 66 by guez, Thu Sep 20 13:00:41 2012 UTC
# Line 1  Line 1 
1  SUBROUTINE integrd(nq,vcovm1,ucovm1,tetam1,psm1,massem1,dv,du,dteta,dq,dp, &  module integrd_m
      vcov,ucov,teta,q,ps,masse,phis,finvmaold,leapf, dt)  
   
   ! From dyn3d/integrd.F,v 1.1.1.1 2004/05/19 12:53:05  
   !   Auteur:  P. Le Van                                                    
   !   objet:                                                                
   !   Incrementation des tendances dynamiques                              
   
   USE dimens_m, ONLY : iim, llm  
   USE paramet_m, ONLY : iip1, iip2, ijp1llm, ip1jm, ip1jmp1, jjp1, llmp1  
   USE comvert, ONLY : ap, bp  
   USE comgeom, ONLY : aire, apoln, apols  
   USE pression_m, ONLY : pression  
   USE filtreg_m, ONLY : filtreg  
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    !   Arguments:                                                            contains
   
   INTEGER nq  
   
   REAL vcov(ip1jm,llm), ucov(ip1jmp1,llm), teta(ip1jmp1,llm)  
   REAL q(ip1jmp1,llm,nq)  
   REAL ps(ip1jmp1), masse(ip1jmp1,llm), phis(ip1jmp1)  
   
   REAL vcovm1(ip1jm,llm), ucovm1(ip1jmp1,llm)  
   REAL tetam1(ip1jmp1,llm), psm1(ip1jmp1), massem1(ip1jmp1,llm)  
   
   REAL dv(ip1jm,llm), du(ip1jmp1,llm)  
   REAL dteta(ip1jmp1,llm), dp(ip1jmp1)  
   REAL dq(ip1jmp1,llm,nq), finvmaold(ip1jmp1,llm)  
   LOGICAL, INTENT (IN) :: leapf  
   real, intent(in):: dt  
   
   !   Local:                                                                
   
   REAL vscr(ip1jm), uscr(ip1jmp1), hscr(ip1jmp1), pscr(ip1jmp1)  
   REAL massescr(ip1jmp1,llm), finvmasse(ip1jmp1,llm)  
   REAL p(ip1jmp1,llmp1)  
   REAL tpn, tps, tppn(iim), tpps(iim)  
   REAL qpn, qps, qppn(iim), qpps(iim)  
   REAL deltap(ip1jmp1,llm)  
   
   INTEGER l, ij, iq  
   
   REAL ssum  
   
   !-----------------------------------------------------------------------  
   
   DO l = 1, llm  
      DO ij = 1, iip1  
         ucov(ij,l) = 0.  
         ucov(ij+ip1jm,l) = 0.  
         uscr(ij) = 0.  
         uscr(ij+ip1jm) = 0.  
      END DO  
   END DO  
   
   
   !    ............    integration  de       ps         ..............      
   
   CALL scopy(ip1jmp1*llm,masse,1,massescr,1)  
   
   DO ij = 1, ip1jmp1  
      pscr(ij) = ps(ij)  
      ps(ij) = psm1(ij) + dt*dp(ij)  
   END DO  
   
   DO ij = 1, ip1jmp1  
      IF (ps(ij)<0.) THEN  
         PRINT *, ' Au point ij = ', ij, ' , pression sol neg. ', ps(ij)  
         STOP 'integrd'  
      END IF  
   END DO  
   
   DO ij = 1, iim  
      tppn(ij) = aire(ij)*ps(ij)  
      tpps(ij) = aire(ij+ip1jm)*ps(ij+ip1jm)  
   END DO  
   tpn = ssum(iim,tppn,1)/apoln  
   tps = ssum(iim,tpps,1)/apols  
   DO ij = 1, iip1  
      ps(ij) = tpn  
      ps(ij+ip1jm) = tps  
   END DO  
   
   !  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 .  
   
   CALL pression(ip1jmp1,ap,bp,ps,p)  
   CALL massdair(p,masse)  
   
   CALL scopy(ijp1llm,masse,1,finvmasse,1)  
   CALL filtreg(finvmasse,jjp1,llm,-2,2,.TRUE.,1)  
   
   
   !    ............   integration  de  ucov, vcov,  h     ..............    
   
   DO  l = 1, llm  
   
      DO ij = iip2, ip1jm  
         uscr(ij) = ucov(ij,l)  
         ucov(ij,l) = ucovm1(ij,l) + dt*du(ij,l)  
      END DO  
   
      DO ij = 1, ip1jm  
         vscr(ij) = vcov(ij,l)  
         vcov(ij,l) = vcovm1(ij,l) + dt*dv(ij,l)  
      END DO  
   
      DO ij = 1, ip1jmp1  
         hscr(ij) = teta(ij,l)  
         teta(ij,l) = tetam1(ij,l)*massem1(ij,l)/masse(ij,l) + &  
              dt*dteta(ij,l)/masse(ij,l)  
      END DO  
   
      !   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    .  
   
   
      DO ij = 1, iim  
         tppn(ij) = aire(ij)*teta(ij,l)  
         tpps(ij) = aire(ij+ip1jm)*teta(ij+ip1jm,l)  
      END DO  
      tpn = ssum(iim,tppn,1)/apoln  
      tps = ssum(iim,tpps,1)/apols  
   
      DO ij = 1, iip1  
         teta(ij,l) = tpn  
         teta(ij+ip1jm,l) = tps  
      END DO  
   
   
      IF (leapf) THEN  
         CALL scopy(ip1jmp1,uscr(1),1,ucovm1(1,l),1)  
         CALL scopy(ip1jm,vscr(1),1,vcovm1(1,l),1)  
         CALL scopy(ip1jmp1,hscr(1),1,tetam1(1,l),1)  
      END IF  
   
   END DO  
   
   DO l = 1, llm  
      DO ij = 1, ip1jmp1  
         deltap(ij,l) = p(ij,l) - p(ij,l+1)  
      END DO  
   END DO  
   
   CALL qminimum(q,nq,deltap)  
   
   !    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....  
   
   
   DO iq = 1, nq  
      DO l = 1, llm  
   
         DO ij = 1, iim  
            qppn(ij) = aire(ij)*q(ij,l,iq)  
            qpps(ij) = aire(ij+ip1jm)*q(ij+ip1jm,l,iq)  
         END DO  
         qpn = ssum(iim,qppn,1)/apoln  
         qps = ssum(iim,qpps,1)/apols  
   
         DO ij = 1, iip1  
            q(ij,l,iq) = qpn  
            q(ij+ip1jm,l,iq) = qps  
         END DO  
   
      END DO  
   END DO  
   
   
   CALL scopy(ijp1llm,finvmasse,1,finvmaold,1)  
6    
7      SUBROUTINE integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, dudyn, &
8           dteta, dp, vcov, ucov, teta, q, ps, masse, finvmaold, dt, leapf)
9    
10    !     .....   FIN  de l'integration  de   q    .......                        ! From dyn3d/integrd.F, version 1.1.1.1 2004/05/19 12:53:05
11        ! Author: P. Le Van
12        ! Objet: incrémentation des tendances dynamiques
13    
14        USE disvert_m, ONLY : ap, bp
15        USE comgeom, ONLY : aire, apoln, apols
16        USE dimens_m, ONLY : iim, jjm, llm
17        USE filtreg_m, ONLY : filtreg
18        use nr_util, only: assert
19        USE paramet_m, ONLY : iip1, iip2, ip1jm, ip1jmp1, jjp1, llmp1
20    
21        ! Arguments:
22    
23        REAL vcov(ip1jm, llm), ucov((iim + 1) * (jjm + 1), llm)
24        real, intent(inout):: teta((iim + 1) * (jjm + 1), llm)
25        REAL q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nq)
26        REAL, intent(inout):: ps((iim + 1) * (jjm + 1))
27        REAL masse((iim + 1) * (jjm + 1), llm)
28    
29        REAL vcovm1(ip1jm, llm), ucovm1((iim + 1) * (jjm + 1), llm)
30        REAL, intent(inout):: tetam1((iim + 1) * (jjm + 1), llm)
31        REAL, intent(inout):: psm1((iim + 1) * (jjm + 1))
32        real massem1((iim + 1) * (jjm + 1), llm)
33    
34        REAL dv(ip1jm, llm), dudyn((iim + 1) * (jjm + 1), llm)
35        REAL dteta((iim + 1) * (jjm + 1), llm), dp((iim + 1) * (jjm + 1))
36        REAL finvmaold((iim + 1) * (jjm + 1), llm)
37        LOGICAL, INTENT (IN) :: leapf
38        real, intent(in):: dt
39    
40        ! Local variables:
41    
42        INTEGER nq
43        REAL vscr(ip1jm), uscr((iim + 1) * (jjm + 1)), hscr((iim + 1) * (jjm + 1))
44        real pscr((iim + 1) * (jjm + 1))
45        REAL massescr((iim + 1) * (jjm + 1), llm)
46        real finvmasse((iim + 1) * (jjm + 1), llm)
47        REAL p((iim + 1) * (jjm + 1), llmp1)
48        REAL tpn, tps, tppn(iim), tpps(iim)
49        REAL qpn, qps, qppn(iim), qpps(iim)
50        REAL deltap((iim + 1) * (jjm + 1), llm)
51    
52        INTEGER l, ij, iq
53    
54        REAL ssum
55    
56        !-----------------------------------------------------------------------
57    
58        call assert(size(q, 1) == iim + 1, size(q, 2) == jjm + 1, &
59             size(q, 3) == llm, "integrd")
60        nq = size(q, 4)
61    
62        DO l = 1, llm
63           DO ij = 1, iip1
64              ucov(ij, l) = 0.
65              ucov(ij+ip1jm, l) = 0.
66              uscr(ij) = 0.
67              uscr(ij+ip1jm) = 0.
68           END DO
69        END DO
70    
71        massescr = masse
72    
73        ! Integration de ps :
74    
75        pscr = ps
76        ps = psm1 + dt * dp
77    
78        DO ij = 1, (iim + 1) * (jjm + 1)
79           IF (ps(ij) < 0.) THEN
80              PRINT *, 'integrd: au point ij = ', ij, &
81                   ', negative surface pressure ', ps(ij)
82              STOP 1
83           END IF
84        END DO
85    
86        DO ij = 1, iim
87           tppn(ij) = aire(ij)*ps(ij)
88           tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
89        END DO
90        tpn = ssum(iim, tppn, 1)/apoln
91        tps = ssum(iim, tpps, 1)/apols
92        DO ij = 1, iip1
93           ps(ij) = tpn
94           ps(ij+ip1jm) = tps
95        END DO
96    
97        ! Calcul de la nouvelle masse d'air au dernier temps integre t+1
98    
99        forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
100        CALL massdair(p, masse)
101    
102        finvmasse = masse
103        CALL filtreg(finvmasse, jjp1, llm, -2, 2, .TRUE.)
104    
105        ! integration de ucov, vcov, h
106    
107        DO l = 1, llm
108           DO ij = iip2, ip1jm
109              uscr(ij) = ucov(ij, l)
110              ucov(ij, l) = ucovm1(ij, l) + dt*dudyn(ij, l)
111           END DO
112    
113           DO ij = 1, ip1jm
114              vscr(ij) = vcov(ij, l)
115              vcov(ij, l) = vcovm1(ij, l) + dt*dv(ij, l)
116           END DO
117    
118           hscr = teta(:, l)
119           teta(:, l) = tetam1(:, l) * massem1(:, l) / masse(:, l) &
120                + dt * dteta(:, l) / masse(:, l)
121    
122           ! Calcul de la valeur moyenne, unique aux poles pour teta
123    
124           DO ij = 1, iim
125              tppn(ij) = aire(ij)*teta(ij, l)
126              tpps(ij) = aire(ij+ip1jm)*teta(ij+ip1jm, l)
127           END DO
128           tpn = ssum(iim, tppn, 1)/apoln
129           tps = ssum(iim, tpps, 1)/apols
130    
131           DO ij = 1, iip1
132              teta(ij, l) = tpn
133              teta(ij+ip1jm, l) = tps
134           END DO
135    
136           IF (leapf) THEN
137              ucovm1(:, l)  =uscr
138              vcovm1(:, l) = vscr
139              tetam1(:, l) = hscr
140           END IF
141        END DO
142    
143        DO l = 1, llm
144           DO ij = 1, (iim + 1) * (jjm + 1)
145              deltap(ij, l) = p(ij, l) - p(ij, l+1)
146           END DO
147        END DO
148    
149        CALL qminimum(q, nq, deltap)
150    
151        ! Calcul de la valeur moyenne, unique aux poles pour q
152    
153        DO iq = 1, nq
154           DO l = 1, llm
155              DO ij = 1, iim
156                 qppn(ij) = aire(ij)*q(ij, 1, l, iq)
157                 qpps(ij) = aire(ij+ip1jm)*q(ij, jjm + 1, l, iq)
158              END DO
159              qpn = ssum(iim, qppn, 1)/apoln
160              qps = ssum(iim, qpps, 1)/apols
161    
162              DO ij = 1, iip1
163                 q(ij, 1, l, iq) = qpn
164                 q(ij, jjm + 1, l, iq) = qps
165              END DO
166           END DO
167        END DO
168    
169        finvmaold = finvmasse
170    
171        ! Fin de l'integration de q
172    
173        IF (leapf) THEN
174           psm1 = pscr
175           massem1 = massescr
176        END IF
177    
178    IF (leapf) THEN    END SUBROUTINE integrd
      CALL scopy(ip1jmp1,pscr,1,psm1,1)  
      CALL scopy(ip1jmp1*llm,massescr,1,massem1,1)  
   END IF  
179    
180  END SUBROUTINE integrd  end module integrd_m
Legend:
Removed from v.28  
changed lines
  Added in v.66
  ViewVC Help
Powered by ViewVC 1.1.21