/[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 260 by guez, Tue Mar 6 17:18:33 2018 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  
   
6    
7    CALL scopy(ijp1llm,finvmasse,1,finvmaold,1)    SUBROUTINE integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, &
8           dp, vcov, ucov, teta, q, ps, masse, dt, leapf)
9    
10        ! 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 massdair_m, only: massdair
18        use nr_util, only: assert
19        USE paramet_m, ONLY : iip1, iip2, ip1jm, llmp1
20        use qminimum_m, only: qminimum
21    
22        REAL vcovm1(ip1jm, llm), ucovm1((iim + 1) * (jjm + 1), llm)
23        REAL, intent(inout):: tetam1(iim + 1, jjm + 1, llm)
24        REAL, intent(inout):: psm1((iim + 1) * (jjm + 1))
25        real, intent(inout):: massem1(iim + 1, jjm + 1, llm)
26        REAL, intent(in):: dv(ip1jm, llm), du((iim + 1) * (jjm + 1), llm)
27        REAL, intent(in):: dteta(iim + 1, jjm + 1, llm), dp((iim + 1) * (jjm + 1))
28        REAL, intent(inout):: vcov(ip1jm, llm), ucov((iim + 1) * (jjm + 1), llm)
29        real, intent(inout):: teta(iim + 1, jjm + 1, llm)
30        REAL q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nq)
31        REAL, intent(inout):: ps((iim + 1) * (jjm + 1))
32        REAL, intent(inout):: masse(iim + 1, jjm + 1, llm)
33        real, intent(in):: dt ! time step, in s
34        LOGICAL, INTENT (IN) :: leapf
35    
36        ! Local:
37        REAL finvmaold(iim + 1, jjm + 1, llm)
38        INTEGER nq
39        REAL vscr(ip1jm), uscr((iim + 1) * (jjm + 1)), hscr(iim + 1, jjm + 1)
40        real pscr((iim + 1) * (jjm + 1))
41        REAL p((iim + 1) * (jjm + 1), llmp1)
42        REAL tpn, tps, tppn(iim), tpps(iim)
43        REAL deltap((iim + 1) * (jjm + 1), llm)
44        INTEGER l, ij, iq
45    
46        !-----------------------------------------------------------------------
47    
48        call assert(size(q, 1) == iim + 1, size(q, 2) == jjm + 1, &
49             size(q, 3) == llm, "integrd")
50        nq = size(q, 4)
51    
52        DO l = 1, llm
53           DO ij = 1, iip1
54              ucov(ij, l) = 0.
55              ucov(ij+ip1jm, l) = 0.
56              uscr(ij) = 0.
57              uscr(ij+ip1jm) = 0.
58           END DO
59        END DO
60    
61        ! Integration de ps :
62    
63        pscr = ps
64        ps = psm1 + dt * dp
65    
66        DO ij = 1, (iim + 1) * (jjm + 1)
67           IF (ps(ij) < 0.) THEN
68              PRINT *, 'integrd: au point ij = ', ij, &
69                   ', negative surface pressure ', ps(ij)
70              STOP 1
71           END IF
72        END DO
73    
74        DO ij = 1, iim
75           tppn(ij) = aire(ij) * ps(ij)
76           tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
77        END DO
78        tpn = sum(tppn)/apoln
79        tps = sum(tpps)/apols
80        DO ij = 1, iip1
81           ps(ij) = tpn
82           ps(ij+ip1jm) = tps
83        END DO
84    
85        ! Calcul de la nouvelle masse d'air au dernier temps integre t+1
86    
87        forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
88        CALL massdair(p, finvmaold)
89    
90        ! integration de ucov, vcov, h
91    
92        DO l = 1, llm
93           DO ij = iip2, ip1jm
94              uscr(ij) = ucov(ij, l)
95              ucov(ij, l) = ucovm1(ij, l) + dt * du(ij, l)
96           END DO
97    
98           DO ij = 1, ip1jm
99              vscr(ij) = vcov(ij, l)
100              vcov(ij, l) = vcovm1(ij, l) + dt * dv(ij, l)
101           END DO
102    
103           hscr = teta(:, :, l)
104           teta(:, :, l) = tetam1(:, :, l) * massem1(:, :, l) / finvmaold(:, :, l) &
105                + dt * dteta(:, :, l) / finvmaold(:, :, l)
106    
107           ! Calcul de la valeur moyenne, unique aux poles pour teta
108           teta(:, 1, l) = sum(aire_2d(:iim, 1) * teta(:iim, 1, l)) / apoln
109           teta(:, jjm + 1, l) = sum(aire_2d(:iim, jjm + 1) &
110                * teta(:iim, jjm + 1, l)) / apols
111    
112           IF (leapf) THEN
113              ucovm1(:, l)  =uscr
114              vcovm1(:, l) = vscr
115              tetam1(:, :, l) = hscr
116           END IF
117        END DO
118    
119        DO l = 1, llm
120           DO ij = 1, (iim + 1) * (jjm + 1)
121              deltap(ij, l) = p(ij, l) - p(ij, l+1)
122           END DO
123        END DO
124    
125        CALL qminimum(q, nq, deltap)
126    
127        ! Calcul de la valeur moyenne, unique aux poles pour q
128        DO iq = 1, nq
129           DO l = 1, llm
130              q(:, 1, l, iq) = sum(aire_2d(:iim, 1) * q(:iim, 1, l, iq)) / apoln
131              q(:, jjm + 1, l, iq) = sum(aire_2d(:iim, jjm + 1) &
132                   * q(:iim, jjm + 1, l, iq)) / apols
133           END DO
134        END DO
135    
136        ! Fin de l'integration de q
137    
138        IF (leapf) THEN
139           psm1 = pscr
140           massem1 = masse
141        END IF
142    
143    !     .....   FIN  de l'integration  de   q    .......                        masse = finvmaold
144    
145    IF (leapf) THEN    END SUBROUTINE integrd
      CALL scopy(ip1jmp1,pscr,1,psm1,1)  
      CALL scopy(ip1jmp1*llm,massescr,1,massem1,1)  
   END IF  
146    
147  END SUBROUTINE integrd  end module integrd_m
Legend:
Removed from v.28  
changed lines
  Added in v.260
  ViewVC Help
Powered by ViewVC 1.1.21