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

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