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

Diff of /trunk/dyn3d/integrd.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/integrd.f90 revision 28 by guez, Fri Mar 26 18:33:04 2010 UTC trunk/dyn3d/integrd.f revision 90 by guez, Wed Mar 12 21:16:36 2014 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, dteta, &
8           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        use qminimum_m, only: qminimum
22    
23        REAL vcovm1(ip1jm, llm), ucovm1((iim + 1) * (jjm + 1), llm)
24        REAL, intent(inout):: tetam1((iim + 1) * (jjm + 1), llm)
25        REAL, intent(inout):: psm1((iim + 1) * (jjm + 1))
26        real massem1((iim + 1) * (jjm + 1), llm)
27        REAL dv(ip1jm, llm), dudyn((iim + 1) * (jjm + 1), llm)
28        REAL dteta((iim + 1) * (jjm + 1), llm), dp((iim + 1) * (jjm + 1))
29        REAL vcov(ip1jm, llm), ucov((iim + 1) * (jjm + 1), llm)
30        real, intent(inout):: teta((iim + 1) * (jjm + 1), llm)
31        REAL q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nq)
32        REAL, intent(inout):: ps((iim + 1) * (jjm + 1))
33        REAL masse((iim + 1) * (jjm + 1), llm)
34        REAL finvmaold((iim + 1) * (jjm + 1), llm)
35        real, intent(in):: dt
36        LOGICAL, INTENT (IN) :: leapf
37    
38        ! Local:
39        INTEGER nq
40        REAL vscr(ip1jm), uscr((iim + 1) * (jjm + 1)), hscr((iim + 1) * (jjm + 1))
41        real pscr((iim + 1) * (jjm + 1))
42        REAL massescr((iim + 1) * (jjm + 1), llm)
43        real finvmasse((iim + 1) * (jjm + 1), llm)
44        REAL p((iim + 1) * (jjm + 1), llmp1)
45        REAL tpn, tps, tppn(iim), tpps(iim)
46        REAL qpn, qps, qppn(iim), qpps(iim)
47        REAL deltap((iim + 1) * (jjm + 1), llm)
48        INTEGER l, ij, iq
49    
50        !-----------------------------------------------------------------------
51    
52        call assert(size(q, 1) == iim + 1, size(q, 2) == jjm + 1, &
53             size(q, 3) == llm, "integrd")
54        nq = size(q, 4)
55    
56        DO l = 1, llm
57           DO ij = 1, iip1
58              ucov(ij, l) = 0.
59              ucov(ij+ip1jm, l) = 0.
60              uscr(ij) = 0.
61              uscr(ij+ip1jm) = 0.
62           END DO
63        END DO
64    
65        massescr = masse
66    
67        ! Integration de ps :
68    
69        pscr = ps
70        ps = psm1 + dt * dp
71    
72        DO ij = 1, (iim + 1) * (jjm + 1)
73           IF (ps(ij) < 0.) THEN
74              PRINT *, 'integrd: au point ij = ', ij, &
75                   ', negative surface pressure ', ps(ij)
76              STOP 1
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 = sum(tppn)/apoln
85        tps = sum(tpps)/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        finvmasse = masse
97        CALL filtreg(finvmasse, jjp1, llm, -2, 2, .TRUE.)
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*dudyn(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           hscr = teta(:, l)
113           teta(:, l) = tetam1(:, l) * massem1(:, l) / masse(:, l) &
114                + dt * dteta(:, l) / masse(:, l)
115    
116           ! Calcul de la valeur moyenne, unique aux poles pour teta
117    
118           DO ij = 1, iim
119              tppn(ij) = aire(ij)*teta(ij, l)
120              tpps(ij) = aire(ij+ip1jm)*teta(ij+ip1jm, l)
121           END DO
122           tpn = sum(tppn)/apoln
123           tps = sum(tpps)/apols
124    
125           DO ij = 1, iip1
126              teta(ij, l) = tpn
127              teta(ij+ip1jm, l) = tps
128           END DO
129    
130           IF (leapf) THEN
131              ucovm1(:, l)  =uscr
132              vcovm1(:, l) = vscr
133              tetam1(:, l) = hscr
134           END IF
135        END DO
136    
137        DO l = 1, llm
138           DO ij = 1, (iim + 1) * (jjm + 1)
139              deltap(ij, l) = p(ij, l) - p(ij, l+1)
140           END DO
141        END DO
142    
143        CALL qminimum(q, nq, deltap)
144    
145        ! Calcul de la valeur moyenne, unique aux poles pour q
146    
147        DO iq = 1, nq
148           DO l = 1, llm
149              DO ij = 1, iim
150                 qppn(ij) = aire(ij)*q(ij, 1, l, iq)
151                 qpps(ij) = aire(ij+ip1jm)*q(ij, jjm + 1, l, iq)
152              END DO
153              qpn = sum(qppn)/apoln
154              qps = sum(qpps)/apols
155    
156              DO ij = 1, iip1
157                 q(ij, 1, l, iq) = qpn
158                 q(ij, jjm + 1, l, iq) = qps
159              END DO
160           END DO
161        END DO
162    
163        finvmaold = finvmasse
164    
165        ! Fin de l'integration de q
166    
167        IF (leapf) THEN
168           psm1 = pscr
169           massem1 = massescr
170        END IF
171    
172    IF (leapf) THEN    END SUBROUTINE integrd
      CALL scopy(ip1jmp1,pscr,1,psm1,1)  
      CALL scopy(ip1jmp1*llm,massescr,1,massem1,1)  
   END IF  
173    
174  END SUBROUTINE integrd  end module integrd_m
Legend:
Removed from v.28  
changed lines
  Added in v.90
  ViewVC Help
Powered by ViewVC 1.1.21