/[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

trunk/libf/dyn3d/integrd.f90 revision 28 by guez, Fri Mar 26 18:33:04 2010 UTC trunk/dyn3d/integrd.f revision 107 by guez, Thu Sep 11 15:09:15 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, aire_2d, 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, intent(in):: dv(ip1jm, llm), dudyn((iim + 1) * (jjm + 1), llm)
28        REAL dteta(iim + 1, jjm + 1, llm), dp((iim + 1) * (jjm + 1))
29        REAL, intent(inout):: 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, intent(inout):: masse(iim + 1, jjm + 1, llm)
34        REAL finvmaold(iim + 1, jjm + 1, llm)
35        real, intent(in):: dt ! time step, in s
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 deltap((iim + 1) * (jjm + 1), llm)
47        INTEGER l, ij, iq
48    
49        !-----------------------------------------------------------------------
50    
51        call assert(size(q, 1) == iim + 1, size(q, 2) == jjm + 1, &
52             size(q, 3) == llm, "integrd")
53        nq = size(q, 4)
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        massescr = masse
65    
66        ! Integration de ps :
67    
68        pscr = ps
69        ps = psm1 + dt * dp
70    
71        DO ij = 1, (iim + 1) * (jjm + 1)
72           IF (ps(ij) < 0.) THEN
73              PRINT *, 'integrd: au point ij = ', ij, &
74                   ', negative surface pressure ', ps(ij)
75              STOP 1
76           END IF
77        END DO
78    
79        DO ij = 1, iim
80           tppn(ij) = aire(ij) * ps(ij)
81           tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
82        END DO
83        tpn = sum(tppn)/apoln
84        tps = sum(tpps)/apols
85        DO ij = 1, iip1
86           ps(ij) = tpn
87           ps(ij+ip1jm) = tps
88        END DO
89    
90        ! Calcul de la nouvelle masse d'air au dernier temps integre t+1
91    
92        forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
93        CALL massdair(p, masse)
94    
95        finvmasse = masse
96        CALL filtreg(finvmasse, direct = .false., intensive = .false.)
97    
98        ! integration de ucov, vcov, h
99    
100        DO l = 1, llm
101           DO ij = iip2, ip1jm
102              uscr(ij) = ucov(ij, l)
103              ucov(ij, l) = ucovm1(ij, l) + dt * dudyn(ij, l)
104           END DO
105    
106           DO ij = 1, ip1jm
107              vscr(ij) = vcov(ij, l)
108              vcov(ij, l) = vcovm1(ij, l) + dt * dv(ij, l)
109           END DO
110    
111           hscr = teta(:, :, l)
112           teta(:, :, l) = tetam1(:, :, l) * massem1(:, :, l) / masse(:, :, l) &
113                + dt * dteta(:, :, l) / masse(:, :, l)
114    
115           ! Calcul de la valeur moyenne, unique aux poles pour teta
116           teta(:, 1, l) = sum(aire_2d(:iim, 1) * teta(:iim, 1, l)) / apoln
117           teta(:, jjm + 1, l) = sum(aire_2d(:iim, jjm + 1) &
118                * teta(:iim, jjm + 1, l)) / apols
119    
120           IF (leapf) THEN
121              ucovm1(:, l)  =uscr
122              vcovm1(:, l) = vscr
123              tetam1(:, :, l) = hscr
124           END IF
125        END DO
126    
127        DO l = 1, llm
128           DO ij = 1, (iim + 1) * (jjm + 1)
129              deltap(ij, l) = p(ij, l) - p(ij, l+1)
130           END DO
131        END DO
132    
133        CALL qminimum(q, nq, deltap)
134    
135        ! Calcul de la valeur moyenne, unique aux poles pour q
136        DO iq = 1, nq
137           DO l = 1, llm
138              q(:, 1, l, iq) = sum(aire_2d(:iim, 1) * q(:iim, 1, l, iq)) / apoln
139              q(:, jjm + 1, l, iq) = sum(aire_2d(:iim, jjm + 1) &
140                   * q(:iim, jjm + 1, l, iq)) / apols
141           END DO
142        END DO
143    
144        finvmaold = finvmasse
145    
146        ! Fin de l'integration de q
147    
148        IF (leapf) THEN
149           psm1 = pscr
150           massem1 = massescr
151        END IF
152    
153    IF (leapf) THEN    END SUBROUTINE integrd
      CALL scopy(ip1jmp1,pscr,1,psm1,1)  
      CALL scopy(ip1jmp1*llm,massescr,1,massem1,1)  
   END IF  
154    
155  END SUBROUTINE integrd  end module integrd_m
Legend:
Removed from v.28  
changed lines
  Added in v.107
  ViewVC Help
Powered by ViewVC 1.1.21