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

Legend:
Removed from v.28  
changed lines
  Added in v.39

  ViewVC Help
Powered by ViewVC 1.1.21