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

Diff of /trunk/dyn3d/integrd.f90

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

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