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

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

  ViewVC Help
Powered by ViewVC 1.1.21