/[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 47 by guez, Fri Jul 1 15:00:48 2011 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 comvert, 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 tetam1((iim + 1) * (jjm + 1), llm), psm1((iim + 1) * (jjm + 1))
31        real massem1((iim + 1) * (jjm + 1), llm)
32    
33        REAL dv(ip1jm, llm), dudyn((iim + 1) * (jjm + 1), llm)
34        REAL dteta((iim + 1) * (jjm + 1), llm), dp((iim + 1) * (jjm + 1))
35        REAL finvmaold((iim + 1) * (jjm + 1), llm)
36        LOGICAL, INTENT (IN) :: leapf
37        real, intent(in):: dt
38    
39        ! Local variables:
40    
41        INTEGER nq
42        REAL vscr(ip1jm), uscr((iim + 1) * (jjm + 1)), hscr((iim + 1) * (jjm + 1))
43        real pscr((iim + 1) * (jjm + 1))
44        REAL massescr((iim + 1) * (jjm + 1), llm)
45        real finvmasse((iim + 1) * (jjm + 1), llm)
46        REAL p((iim + 1) * (jjm + 1), llmp1)
47        REAL tpn, tps, tppn(iim), tpps(iim)
48        REAL qpn, qps, qppn(iim), qpps(iim)
49        REAL deltap((iim + 1) * (jjm + 1), llm)
50    
51        INTEGER l, ij, iq
52    
53        REAL ssum
54    
55        !-----------------------------------------------------------------------
56    
57        call assert(size(q, 1) == iim + 1, size(q, 2) == jjm + 1, &
58             size(q, 3) == llm, "integrd")
59        nq = size(q, 4)
60    
61        DO l = 1, llm
62           DO ij = 1, iip1
63              ucov(ij, l) = 0.
64              ucov(ij+ip1jm, l) = 0.
65              uscr(ij) = 0.
66              uscr(ij+ip1jm) = 0.
67           END DO
68        END DO
69    
70        massescr = masse
71    
72        ! Integration de ps :
73    
74        pscr = ps
75        ps = psm1 + dt * dp
76    
77        DO ij = 1, (iim + 1) * (jjm + 1)
78           IF (ps(ij) < 0.) THEN
79              PRINT *, 'integrd: au point ij = ', ij, &
80                   ', negative surface pressure ', ps(ij)
81              STOP 1
82           END IF
83        END DO
84    
85        DO ij = 1, iim
86           tppn(ij) = aire(ij)*ps(ij)
87           tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
88        END DO
89        tpn = ssum(iim, tppn, 1)/apoln
90        tps = ssum(iim, tpps, 1)/apols
91        DO ij = 1, iip1
92           ps(ij) = tpn
93           ps(ij+ip1jm) = tps
94        END DO
95    
96        ! Calcul de la nouvelle masse d'air au dernier temps integre t+1
97    
98        forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
99        CALL massdair(p, masse)
100    
101        finvmasse = masse
102        CALL filtreg(finvmasse, jjp1, llm, -2, 2, .TRUE., 1)
103    
104        ! integration de ucov, vcov, h
105    
106        DO l = 1, llm
107           DO ij = iip2, ip1jm
108              uscr(ij) = ucov(ij, l)
109              ucov(ij, l) = ucovm1(ij, l) + dt*dudyn(ij, l)
110           END DO
111    
112           DO ij = 1, ip1jm
113              vscr(ij) = vcov(ij, l)
114              vcov(ij, l) = vcovm1(ij, l) + dt*dv(ij, l)
115           END DO
116    
117           DO ij = 1, (iim + 1) * (jjm + 1)
118              hscr(ij) = teta(ij, l)
119              teta(ij, l) = tetam1(ij, l) * massem1(ij, l) / masse(ij, l) &
120                   + dt * dteta(ij, l) / masse(ij, l)
121           END DO
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              CALL scopy((iim + 1) * (jjm + 1), uscr(1), 1, ucovm1(1, l), 1)
139              CALL scopy(ip1jm, vscr(1), 1, vcovm1(1, l), 1)
140              CALL scopy((iim + 1) * (jjm + 1), hscr(1), 1, tetam1(1, l), 1)
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        CALL scopy((iim + 1) * (jjm + 1) * llm, finvmasse, 1, finvmaold, 1)
171    
172        ! Fin de l'integration de q
173    
174        IF (leapf) THEN
175           CALL scopy((iim + 1) * (jjm + 1), pscr, 1, psm1, 1)
176           CALL scopy((iim + 1) * (jjm + 1)*llm, massescr, 1, massem1, 1)
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.47

  ViewVC Help
Powered by ViewVC 1.1.21