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

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

  ViewVC Help
Powered by ViewVC 1.1.21