/[lmdze]/trunk/Sources/dyn3d/pentes_ini.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/pentes_ini.f

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

trunk/libf/dyn3d/pentes_ini.f revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/dyn3d/pentes_ini.f revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC
# Line 1  Line 1 
 !  
 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/pentes_ini.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $  
 !  
       SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)  
       use dimens_m  
       use paramet_m  
       use comconst  
       use comvert  
       use comgeom  
       IMPLICIT NONE  
   
 c=======================================================================  
 c   Adaptation LMDZ:  A.Armengaud (LGGE)  
 c   ----------------  
 c  
 c   ********************************************************************  
 c   Transport des traceurs par la methode des pentes  
 c   ********************************************************************  
 c   Reference possible : Russel. G.L., Lerner J.A.:  
 c         A new Finite-Differencing Scheme for Traceur Transport  
 c         Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81  
 c   ********************************************************************  
 c   q,w,masse,pbaru et pbarv  
 c                      sont des arguments d'entree  pour le s-pg ....  
 c  
 c=======================================================================  
   
   
   
 c   Arguments:  
 c   ----------  
       integer mode  
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )  
       REAL q( iip1,jjp1,llm,0:3)  
       REAL w( ip1jmp1,llm )  
       REAL masse( iip1,jjp1,llm)  
 c   Local:  
 c   ------  
       LOGICAL limit  
       REAL sm ( iip1,jjp1, llm )  
       REAL s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )  
       REAL sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )  
       real masn,mass,zz  
       INTEGER i,j,l,iq  
   
 c  modif Fred 24 03 96  
   
       real sinlon(iip1),sinlondlon(iip1)  
       real coslon(iip1),coslondlon(iip1)  
       save sinlon,coslon,sinlondlon,coslondlon  
       real dyn1,dyn2,dys1,dys2  
       real qpn,qps,dqzpn,dqzps  
       real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)  
       real qmin,zq,pente_max  
 c  
       REAL      SSUM  
       integer ismax,ismin,lati,latf  
       EXTERNAL  SSUM, convflu,ismin,ismax  
       logical first  
       save first  
 c   fin modif  
   
 c      EXTERNAL masskg  
       EXTERNAL advx  
       EXTERNAL advy  
       EXTERNAL advz  
   
 c  modif Fred 24 03 96  
       data first/.true./  
   
       limit = .TRUE.  
       pente_max=2  
 c     if (mode.eq.1.or.mode.eq.3) then  
 c     if (mode.eq.1) then  
       if (mode.ge.1) then  
         lati=2  
         latf=jjm  
       else  
         lati=1  
         latf=jjp1  
       endif  
   
       qmin=0.4995  
       qmin=0.  
       if(first) then  
          print*,'SCHEMA AMONT NOUVEAU'  
          first=.false.  
          do i=2,iip1  
             coslon(i)=cos(rlonv(i))  
             sinlon(i)=sin(rlonv(i))  
             coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi  
             sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi  
             print*,coslondlon(i),sinlondlon(i)  
          enddo  
          coslon(1)=coslon(iip1)  
          coslondlon(1)=coslondlon(iip1)  
          sinlon(1)=sinlon(iip1)  
          sinlondlon(1)=sinlondlon(iip1)  
          print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1)  
          print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1)  
         DO l = 1,llm  
         DO j = 1,jjp1  
          DO i = 1,iip1  
          q ( i,j,l,1 )=0.  
          q ( i,j,l,2 )=0.  
          q ( i,j,l,3 )=0.    
          ENDDO  
          ENDDO  
         ENDDO  
           
       endif  
 c   Fin modif Fred  
   
 c *** q contient les qqtes de traceur avant l'advection  
   
 c *** Affectation des tableaux S a partir de Q  
 c *** Rem : utilisation de SCOPY ulterieurement  
   
        DO l = 1,llm  
         DO j = 1,jjp1  
          DO i = 1,iip1  
              s0( i,j,llm+1-l ) = q ( i,j,l,0 )  
              sx( i,j,llm+1-l ) = q ( i,j,l,1 )  
              sy( i,j,llm+1-l ) = q ( i,j,l,2 )  
              sz( i,j,llm+1-l ) = q ( i,j,l,3 )  
          ENDDO  
         ENDDO  
        ENDDO  
   
 c      PRINT*,'----- S0 just before conversion -------'  
 c      PRINT*,'S0(16,12,1)=',s0(16,12,1)  
 c      PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)  
   
 c *** On calcule la masse d'air en kg  
   
        DO  l = 1,llm  
          DO  j = 1,jjp1  
            DO  i = 1,iip1  
             sm ( i,j,llm+1-l)=masse( i,j,l )  
           ENDDO  
          ENDDO  
        ENDDO  
   
 c *** On converti les champs S en atome (resp. kg)  
 c *** Les routines d'advection traitent les champs  
 c *** a advecter si ces derniers sont en atome (resp. kg)  
 c *** A optimiser !!!  
   
        DO  l = 1,llm  
          DO  j = 1,jjp1  
            DO  i = 1,iip1  
                s0(i,j,l) = s0(i,j,l) * sm ( i,j,l )  
                sx(i,j,l) = sx(i,j,l) * sm ( i,j,l )  
                sy(i,j,l) = sy(i,j,l) * sm ( i,j,l )  
                sz(i,j,l) = sz(i,j,l) * sm ( i,j,l )  
            ENDDO  
          ENDDO  
        ENDDO  
   
 c       ss0 = 0.  
 c       DO l = 1,llm  
 c        DO j = 1,jjp1  
 c         DO i = 1,iim  
 c            ss0 = ss0 + s0 ( i,j,l )  
 c         ENDDO  
 c        ENDDO  
 c       ENDDO  
 c       PRINT*, 'valeur tot s0 avant advection=',ss0  
   
 c *** Appel des subroutines d'advection en X, en Y et en Z  
 c *** Advection avec "time-splitting"  
         
 c-----------------------------------------------------------  
 c      PRINT*,'----- S0 just before ADVX -------'  
 c      PRINT*,'S0(16,12,1)=',s0(16,12,1)  
   
 c-----------------------------------------------------------  
 c      do l=1,llm  
 c         do j=1,jjp1  
 c          do i=1,iip1  
 c             zq=s0(i,j,l)/sm(i,j,l)  
 c            if(zq.lt.qmin)  
 c    ,       print*,'avant advx1, s0(',i,',',j,',',l,')=',zq  
 c          enddo  
 c         enddo  
 c      enddo  
 CCC  
        if(mode.eq.2) then  
           do l=1,llm  
             s0s=0.  
             s0n=0.  
             dyn1=0.  
             dys1=0.  
             dyn2=0.  
             dys2=0.  
             smn=0.  
             sms=0.  
             do i=1,iim  
                smn=smn+sm(i,1,l)  
                sms=sms+sm(i,jjp1,l)  
                s0n=s0n+s0(i,1,l)  
                s0s=s0s+s0(i,jjp1,l)  
                zz=sy(i,1,l)/sm(i,1,l)  
                dyn1=dyn1+sinlondlon(i)*zz  
                dyn2=dyn2+coslondlon(i)*zz  
                zz=sy(i,jjp1,l)/sm(i,jjp1,l)  
                dys1=dys1+sinlondlon(i)*zz  
                dys2=dys2+coslondlon(i)*zz  
             enddo  
             do i=1,iim  
                sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i)  
                sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i)  
             enddo  
             do i=1,iim  
                s0(i,1,l)=s0n/smn+sy(i,1,l)  
                s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)  
             enddo  
   
             s0(iip1,1,l)=s0(1,1,l)  
             s0(iip1,jjp1,l)=s0(1,jjp1,l)  
   
             do i=1,iim  
                sxn(i)=s0(i+1,1,l)-s0(i,1,l)  
                sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)  
 c   on rerentre les masses  
             enddo  
             do i=1,iim  
                sy(i,1,l)=sy(i,1,l)*sm(i,1,l)  
                sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)  
                s0(i,1,l)=s0(i,1,l)*sm(i,1,l)  
                s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)  
             enddo  
             sxn(iip1)=sxn(1)  
             sxs(iip1)=sxs(1)  
             do i=1,iim  
                sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l)  
                sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l)  
             enddo  
             s0(iip1,1,l)=s0(1,1,l)  
             s0(iip1,jjp1,l)=s0(1,jjp1,l)  
             sy(iip1,1,l)=sy(1,1,l)  
             sy(iip1,jjp1,l)=sy(1,jjp1,l)  
             sx(1,1,l)=sx(iip1,1,l)  
             sx(1,jjp1,l)=sx(iip1,jjp1,l)  
           enddo  
       endif  
   
       if (mode.eq.4) then  
          do l=1,llm  
             do i=1,iip1  
                sx(i,1,l)=0.  
                sx(i,jjp1,l)=0.  
                sy(i,1,l)=0.  
                sy(i,jjp1,l)=0.  
             enddo  
          enddo  
       endif  
       call limx(s0,sx,sm,pente_max)  
 c     call minmaxq(zq,1.e33,-1.e33,'avant advx     ')  
        call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)  
 c     call minmaxq(zq,1.e33,-1.e33,'avant advy     ')  
       if (mode.eq.4) then  
          do l=1,llm  
             do i=1,iip1  
                sx(i,1,l)=0.  
                sx(i,jjp1,l)=0.  
                sy(i,1,l)=0.  
                sy(i,jjp1,l)=0.  
             enddo  
          enddo  
       endif  
        call   limy(s0,sy,sm,pente_max)  
        call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )  
 c     call minmaxq(zq,1.e33,-1.e33,'avant advz     ')  
        do j=1,jjp1  
           do i=1,iip1  
              sz(i,j,1)=0.  
              sz(i,j,llm)=0.  
           enddo  
        enddo  
        call limz(s0,sz,sm,pente_max)  
        call advz( limit,dtvr,w,sm,s0,sx,sy,sz )  
       if (mode.eq.4) then  
          do l=1,llm  
             do i=1,iip1  
                sx(i,1,l)=0.  
                sx(i,jjp1,l)=0.  
                sy(i,1,l)=0.  
                sy(i,jjp1,l)=0.  
             enddo  
          enddo  
       endif  
         call limy(s0,sy,sm,pente_max)  
        call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )  
        do l=1,llm  
           do j=1,jjp1  
              sm(iip1,j,l)=sm(1,j,l)  
              s0(iip1,j,l)=s0(1,j,l)  
              sx(iip1,j,l)=sx(1,j,l)  
              sy(iip1,j,l)=sy(1,j,l)  
              sz(iip1,j,l)=sz(1,j,l)  
           enddo  
        enddo  
   
   
 c     call minmaxq(zq,1.e33,-1.e33,'avant advx     ')  
       if (mode.eq.4) then  
          do l=1,llm  
             do i=1,iip1  
                sx(i,1,l)=0.  
                sx(i,jjp1,l)=0.  
                sy(i,1,l)=0.  
                sy(i,jjp1,l)=0.  
             enddo  
          enddo  
       endif  
        call limx(s0,sx,sm,pente_max)  
        call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)  
 c     call minmaxq(zq,1.e33,-1.e33,'apres advx     ')  
 c      do l=1,llm  
 c         do j=1,jjp1  
 c          do i=1,iip1  
 c             zq=s0(i,j,l)/sm(i,j,l)  
 c            if(zq.lt.qmin)  
 c    ,       print*,'apres advx2, s0(',i,',',j,',',l,')=',zq  
 c          enddo  
 c         enddo  
 c      enddo  
 c ***   On repasse les S dans la variable q directement 14/10/94  
 c   On revient a des rapports de melange en divisant par la masse  
   
 c En dehors des poles:  
   
        DO  l = 1,llm  
         DO  j = 1,jjp1  
          DO  i = 1,iim  
              q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)  
              q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)  
              q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)  
              q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)  
          ENDDO  
         ENDDO  
       ENDDO  
   
 c Traitements specifiques au pole  
   
       if(mode.ge.1) then  
       DO l=1,llm  
 c   filtrages aux poles  
          masn=ssum(iim,sm(1,1,l),1)  
          mass=ssum(iim,sm(1,jjp1,l),1)  
          qpn=ssum(iim,s0(1,1,l),1)/masn  
          qps=ssum(iim,s0(1,jjp1,l),1)/mass  
          dqzpn=ssum(iim,sz(1,1,l),1)/masn  
          dqzps=ssum(iim,sz(1,jjp1,l),1)/mass  
          do i=1,iip1  
             q( i,1,llm+1-l,3)=dqzpn  
             q( i,jjp1,llm+1-l,3)=dqzps  
             q( i,1,llm+1-l,0)=qpn  
             q( i,jjp1,llm+1-l,0)=qps  
          enddo  
          if(mode.eq.3) then  
             dyn1=0.  
             dys1=0.  
             dyn2=0.  
             dys2=0.  
             do i=1,iim  
                dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)  
                dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)  
                dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)  
                dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)  
             enddo  
             do i=1,iim  
                q(i,1,llm+1-l,2)=  
      s          (sinlon(i)*dyn1+coslon(i)*dyn2)  
                q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)  
                q(i,jjp1,llm+1-l,2)=  
      s          (sinlon(i)*dys1+coslon(i)*dys2)  
                q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)  
      s         -q(i,jjp1,llm+1-l,2)  
             enddo  
          endif  
          if(mode.eq.1) then  
 c   on filtre les valeurs au bord de la "grande maille pole"  
             dyn1=0.  
             dys1=0.  
             dyn2=0.  
             dys2=0.  
             do i=1,iim  
                zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)  
                dyn1=dyn1+sinlondlon(i)*zz  
                dyn2=dyn2+coslondlon(i)*zz  
                zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)  
                dys1=dys1+sinlondlon(i)*zz  
                dys2=dys2+coslondlon(i)*zz  
             enddo  
             do i=1,iim  
                q(i,1,llm+1-l,2)=  
      s          (sinlon(i)*dyn1+coslon(i)*dyn2)/2.  
                q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)  
                q(i,jjp1,llm+1-l,2)=  
      s          (sinlon(i)*dys1+coslon(i)*dys2)/2.  
                q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)  
      s         -q(i,jjp1,llm+1-l,2)  
             enddo  
             q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)  
             q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)  
   
             do i=1,iim  
                sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)  
                sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)  
             enddo  
             sxn(iip1)=sxn(1)  
             sxs(iip1)=sxs(1)  
             do i=1,iim  
                q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))  
                q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))  
             enddo  
             q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)  
             q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)  
   
          endif  
   
        ENDDO  
        endif  
   
 c bouclage en longitude  
       do iq=0,3  
          do l=1,llm  
             do j=1,jjp1  
                q(iip1,j,l,iq)=q(1,j,l,iq)  
             enddo  
          enddo  
       enddo  
   
 c       PRINT*, ' SORTIE DE PENTES ---  ca peut glisser ....'  
   
         DO l = 1,llm  
          DO j = 1,jjp1  
           DO i = 1,iip1  
                 IF (q(i,j,l,0).lt.0.)  THEN  
 c                    PRINT*,'------------ BIP-----------'  
 c                    PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)  
 c                    PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)  
 c                    PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)  
 c                    PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)  
 c                            PRINT*,' PBL EN SORTIE DE PENTES'  
                      q(i,j,l,0)=0.  
 c                    STOP  
                  ENDIF  
           ENDDO  
          ENDDO  
         ENDDO  
   
 c       PRINT*, '-------------------------------------------'  
           
        do l=1,llm  
           do j=1,jjp1  
            do i=1,iip1  
              if(q(i,j,l,0).lt.qmin)  
      ,       print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)  
            enddo  
           enddo  
        enddo  
       RETURN  
       END  
   
   
   
   
   
   
   
   
   
   
1    
2    ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/pentes_ini.F,v 1.1.1.1 2004/05/19
3    ! 12:53:07 lmdzadmin Exp $
4    
5    SUBROUTINE pentes_ini(q, w, masse, pbaru, pbarv, mode)
6      USE dimens_m
7      USE paramet_m
8      USE comconst
9      USE disvert_m
10      USE comgeom
11      USE nr_util, ONLY: pi
12      IMPLICIT NONE
13    
14      ! =======================================================================
15      ! Adaptation LMDZ:  A.Armengaud (LGGE)
16      ! ----------------
17    
18      ! ********************************************************************
19      ! Transport des traceurs par la methode des pentes
20      ! ********************************************************************
21      ! Reference possible : Russel. G.L., Lerner J.A.:
22      ! A new Finite-Differencing Scheme for Traceur Transport
23      ! Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
24      ! ********************************************************************
25      ! q,w,masse,pbaru et pbarv
26      ! sont des arguments d'entree  pour le s-pg ....
27    
28      ! =======================================================================
29    
30    
31    
32      ! Arguments:
33      ! ----------
34      INTEGER mode
35      REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
36      REAL q(iip1, jjp1, llm, 0:3)
37      REAL w(ip1jmp1, llm)
38      REAL masse(iip1, jjp1, llm)
39      ! Local:
40      ! ------
41      LOGICAL limit
42      REAL sm(iip1, jjp1, llm)
43      REAL s0(iip1, jjp1, llm), sx(iip1, jjp1, llm)
44      REAL sy(iip1, jjp1, llm), sz(iip1, jjp1, llm)
45      REAL masn, mass, zz
46      INTEGER i, j, l, iq
47    
48      ! modif Fred 24 03 96
49    
50      REAL sinlon(iip1), sinlondlon(iip1)
51      REAL coslon(iip1), coslondlon(iip1)
52      SAVE sinlon, coslon, sinlondlon, coslondlon
53      REAL dyn1, dyn2, dys1, dys2
54      REAL qpn, qps, dqzpn, dqzps
55      REAL smn, sms, s0n, s0s, sxn(iip1), sxs(iip1)
56      REAL qmin, zq, pente_max
57    
58      REAL ssum
59      INTEGER ismax, ismin, lati, latf
60      EXTERNAL ssum, convflu, ismin, ismax
61      LOGICAL first
62      SAVE first
63      ! fin modif
64    
65      ! EXTERNAL masskg
66      EXTERNAL advx
67      EXTERNAL advy
68      EXTERNAL advz
69    
70      ! modif Fred 24 03 96
71      DATA first/.TRUE./
72    
73      limit = .TRUE.
74      pente_max = 2
75      ! if (mode.eq.1.or.mode.eq.3) then
76      ! if (mode.eq.1) then
77      IF (mode>=1) THEN
78        lati = 2
79        latf = jjm
80      ELSE
81        lati = 1
82        latf = jjp1
83      END IF
84    
85      qmin = 0.4995
86      qmin = 0.
87      IF (first) THEN
88        PRINT *, 'SCHEMA AMONT NOUVEAU'
89        first = .FALSE.
90        DO i = 2, iip1
91          coslon(i) = cos(rlonv(i))
92          sinlon(i) = sin(rlonv(i))
93          coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
94          sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
95          PRINT *, coslondlon(i), sinlondlon(i)
96        END DO
97        coslon(1) = coslon(iip1)
98        coslondlon(1) = coslondlon(iip1)
99        sinlon(1) = sinlon(iip1)
100        sinlondlon(1) = sinlondlon(iip1)
101        PRINT *, 'sum sinlondlon ', ssum(iim, sinlondlon, 1)/sinlondlon(1)
102        PRINT *, 'sum coslondlon ', ssum(iim, coslondlon, 1)/coslondlon(1)
103        DO l = 1, llm
104          DO j = 1, jjp1
105            DO i = 1, iip1
106              q(i, j, l, 1) = 0.
107              q(i, j, l, 2) = 0.
108              q(i, j, l, 3) = 0.
109            END DO
110          END DO
111        END DO
112    
113      END IF
114    
115      ! *** q contient les qqtes de traceur avant l'advection
116    
117      ! *** Affectation des tableaux S a partir de Q
118      ! *** Rem : utilisation de SCOPY ulterieurement
119    
120      DO l = 1, llm
121        DO j = 1, jjp1
122          DO i = 1, iip1
123            s0(i, j, llm+1-l) = q(i, j, l, 0)
124            sx(i, j, llm+1-l) = q(i, j, l, 1)
125            sy(i, j, llm+1-l) = q(i, j, l, 2)
126            sz(i, j, llm+1-l) = q(i, j, l, 3)
127          END DO
128        END DO
129      END DO
130    
131      ! *** On calcule la masse d'air en kg
132    
133      DO l = 1, llm
134        DO j = 1, jjp1
135          DO i = 1, iip1
136            sm(i, j, llm+1-l) = masse(i, j, l)
137          END DO
138        END DO
139      END DO
140    
141      ! *** On converti les champs S en atome (resp. kg)
142      ! *** Les routines d'advection traitent les champs
143      ! *** a advecter si ces derniers sont en atome (resp. kg)
144      ! *** A optimiser !!!
145    
146      DO l = 1, llm
147        DO j = 1, jjp1
148          DO i = 1, iip1
149            s0(i, j, l) = s0(i, j, l)*sm(i, j, l)
150            sx(i, j, l) = sx(i, j, l)*sm(i, j, l)
151            sy(i, j, l) = sy(i, j, l)*sm(i, j, l)
152            sz(i, j, l) = sz(i, j, l)*sm(i, j, l)
153          END DO
154        END DO
155      END DO
156    
157      ! *** Appel des subroutines d'advection en X, en Y et en Z
158      ! *** Advection avec "time-splitting"
159    
160      IF (mode==2) THEN
161        DO l = 1, llm
162          s0s = 0.
163          s0n = 0.
164          dyn1 = 0.
165          dys1 = 0.
166          dyn2 = 0.
167          dys2 = 0.
168          smn = 0.
169          sms = 0.
170          DO i = 1, iim
171            smn = smn + sm(i, 1, l)
172            sms = sms + sm(i, jjp1, l)
173            s0n = s0n + s0(i, 1, l)
174            s0s = s0s + s0(i, jjp1, l)
175            zz = sy(i, 1, l)/sm(i, 1, l)
176            dyn1 = dyn1 + sinlondlon(i)*zz
177            dyn2 = dyn2 + coslondlon(i)*zz
178            zz = sy(i, jjp1, l)/sm(i, jjp1, l)
179            dys1 = dys1 + sinlondlon(i)*zz
180            dys2 = dys2 + coslondlon(i)*zz
181          END DO
182          DO i = 1, iim
183            sy(i, 1, l) = dyn1*sinlon(i) + dyn2*coslon(i)
184            sy(i, jjp1, l) = dys1*sinlon(i) + dys2*coslon(i)
185          END DO
186          DO i = 1, iim
187            s0(i, 1, l) = s0n/smn + sy(i, 1, l)
188            s0(i, jjp1, l) = s0s/sms - sy(i, jjp1, l)
189          END DO
190    
191          s0(iip1, 1, l) = s0(1, 1, l)
192          s0(iip1, jjp1, l) = s0(1, jjp1, l)
193    
194          DO i = 1, iim
195            sxn(i) = s0(i+1, 1, l) - s0(i, 1, l)
196            sxs(i) = s0(i+1, jjp1, l) - s0(i, jjp1, l)
197            ! on rerentre les masses
198          END DO
199          DO i = 1, iim
200            sy(i, 1, l) = sy(i, 1, l)*sm(i, 1, l)
201            sy(i, jjp1, l) = sy(i, jjp1, l)*sm(i, jjp1, l)
202            s0(i, 1, l) = s0(i, 1, l)*sm(i, 1, l)
203            s0(i, jjp1, l) = s0(i, jjp1, l)*sm(i, jjp1, l)
204          END DO
205          sxn(iip1) = sxn(1)
206          sxs(iip1) = sxs(1)
207          DO i = 1, iim
208            sx(i+1, 1, l) = 0.25*(sxn(i)+sxn(i+1))*sm(i+1, 1, l)
209            sx(i+1, jjp1, l) = 0.25*(sxs(i)+sxs(i+1))*sm(i+1, jjp1, l)
210          END DO
211          s0(iip1, 1, l) = s0(1, 1, l)
212          s0(iip1, jjp1, l) = s0(1, jjp1, l)
213          sy(iip1, 1, l) = sy(1, 1, l)
214          sy(iip1, jjp1, l) = sy(1, jjp1, l)
215          sx(1, 1, l) = sx(iip1, 1, l)
216          sx(1, jjp1, l) = sx(iip1, jjp1, l)
217        END DO
218      END IF
219    
220      IF (mode==4) THEN
221        DO l = 1, llm
222          DO i = 1, iip1
223            sx(i, 1, l) = 0.
224            sx(i, jjp1, l) = 0.
225            sy(i, 1, l) = 0.
226            sy(i, jjp1, l) = 0.
227          END DO
228        END DO
229      END IF
230      CALL limx(s0, sx, sm, pente_max)
231      CALL advx(limit, .5*dtvr, pbaru, sm, s0, sx, sy, sz, lati, latf)
232      IF (mode==4) THEN
233        DO l = 1, llm
234          DO i = 1, iip1
235            sx(i, 1, l) = 0.
236            sx(i, jjp1, l) = 0.
237            sy(i, 1, l) = 0.
238            sy(i, jjp1, l) = 0.
239          END DO
240        END DO
241      END IF
242      CALL limy(s0, sy, sm, pente_max)
243      CALL advy(limit, .5*dtvr, pbarv, sm, s0, sx, sy, sz)
244      DO j = 1, jjp1
245        DO i = 1, iip1
246          sz(i, j, 1) = 0.
247          sz(i, j, llm) = 0.
248        END DO
249      END DO
250      CALL limz(s0, sz, sm, pente_max)
251      CALL advz(limit, dtvr, w, sm, s0, sx, sy, sz)
252      IF (mode==4) THEN
253        DO l = 1, llm
254          DO i = 1, iip1
255            sx(i, 1, l) = 0.
256            sx(i, jjp1, l) = 0.
257            sy(i, 1, l) = 0.
258            sy(i, jjp1, l) = 0.
259          END DO
260        END DO
261      END IF
262      CALL limy(s0, sy, sm, pente_max)
263      CALL advy(limit, .5*dtvr, pbarv, sm, s0, sx, sy, sz)
264      DO l = 1, llm
265        DO j = 1, jjp1
266          sm(iip1, j, l) = sm(1, j, l)
267          s0(iip1, j, l) = s0(1, j, l)
268          sx(iip1, j, l) = sx(1, j, l)
269          sy(iip1, j, l) = sy(1, j, l)
270          sz(iip1, j, l) = sz(1, j, l)
271        END DO
272      END DO
273    
274    
275      IF (mode==4) THEN
276        DO l = 1, llm
277          DO i = 1, iip1
278            sx(i, 1, l) = 0.
279            sx(i, jjp1, l) = 0.
280            sy(i, 1, l) = 0.
281            sy(i, jjp1, l) = 0.
282          END DO
283        END DO
284      END IF
285      CALL limx(s0, sx, sm, pente_max)
286      CALL advx(limit, .5*dtvr, pbaru, sm, s0, sx, sy, sz, lati, latf)
287      ! ***   On repasse les S dans la variable q directement 14/10/94
288      ! On revient a des rapports de melange en divisant par la masse
289    
290      ! En dehors des poles:
291    
292      DO l = 1, llm
293        DO j = 1, jjp1
294          DO i = 1, iim
295            q(i, j, llm+1-l, 0) = s0(i, j, l)/sm(i, j, l)
296            q(i, j, llm+1-l, 1) = sx(i, j, l)/sm(i, j, l)
297            q(i, j, llm+1-l, 2) = sy(i, j, l)/sm(i, j, l)
298            q(i, j, llm+1-l, 3) = sz(i, j, l)/sm(i, j, l)
299          END DO
300        END DO
301      END DO
302    
303      ! Traitements specifiques au pole
304    
305      IF (mode>=1) THEN
306        DO l = 1, llm
307          ! filtrages aux poles
308          masn = ssum(iim, sm(1,1,l), 1)
309          mass = ssum(iim, sm(1,jjp1,l), 1)
310          qpn = ssum(iim, s0(1,1,l), 1)/masn
311          qps = ssum(iim, s0(1,jjp1,l), 1)/mass
312          dqzpn = ssum(iim, sz(1,1,l), 1)/masn
313          dqzps = ssum(iim, sz(1,jjp1,l), 1)/mass
314          DO i = 1, iip1
315            q(i, 1, llm+1-l, 3) = dqzpn
316            q(i, jjp1, llm+1-l, 3) = dqzps
317            q(i, 1, llm+1-l, 0) = qpn
318            q(i, jjp1, llm+1-l, 0) = qps
319          END DO
320          IF (mode==3) THEN
321            dyn1 = 0.
322            dys1 = 0.
323            dyn2 = 0.
324            dys2 = 0.
325            DO i = 1, iim
326              dyn1 = dyn1 + sinlondlon(i)*sy(i, 1, l)/sm(i, 1, l)
327              dyn2 = dyn2 + coslondlon(i)*sy(i, 1, l)/sm(i, 1, l)
328              dys1 = dys1 + sinlondlon(i)*sy(i, jjp1, l)/sm(i, jjp1, l)
329              dys2 = dys2 + coslondlon(i)*sy(i, jjp1, l)/sm(i, jjp1, l)
330            END DO
331            DO i = 1, iim
332              q(i, 1, llm+1-l, 2) = (sinlon(i)*dyn1+coslon(i)*dyn2)
333              q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
334              q(i, jjp1, llm+1-l, 2) = (sinlon(i)*dys1+coslon(i)*dys2)
335              q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
336                q(i, jjp1, llm+1-l, 2)
337            END DO
338          END IF
339          IF (mode==1) THEN
340            ! on filtre les valeurs au bord de la "grande maille pole"
341            dyn1 = 0.
342            dys1 = 0.
343            dyn2 = 0.
344            dys2 = 0.
345            DO i = 1, iim
346              zz = s0(i, 2, l)/sm(i, 2, l) - q(i, 1, llm+1-l, 0)
347              dyn1 = dyn1 + sinlondlon(i)*zz
348              dyn2 = dyn2 + coslondlon(i)*zz
349              zz = q(i, jjp1, llm+1-l, 0) - s0(i, jjm, l)/sm(i, jjm, l)
350              dys1 = dys1 + sinlondlon(i)*zz
351              dys2 = dys2 + coslondlon(i)*zz
352            END DO
353            DO i = 1, iim
354              q(i, 1, llm+1-l, 2) = (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
355              q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
356              q(i, jjp1, llm+1-l, 2) = (sinlon(i)*dys1+coslon(i)*dys2)/2.
357              q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
358                q(i, jjp1, llm+1-l, 2)
359            END DO
360            q(iip1, 1, llm+1-l, 0) = q(1, 1, llm+1-l, 0)
361            q(iip1, jjp1, llm+1-l, 0) = q(1, jjp1, llm+1-l, 0)
362    
363            DO i = 1, iim
364              sxn(i) = q(i+1, 1, llm+1-l, 0) - q(i, 1, llm+1-l, 0)
365              sxs(i) = q(i+1, jjp1, llm+1-l, 0) - q(i, jjp1, llm+1-l, 0)
366            END DO
367            sxn(iip1) = sxn(1)
368            sxs(iip1) = sxs(1)
369            DO i = 1, iim
370              q(i+1, 1, llm+1-l, 1) = 0.25*(sxn(i)+sxn(i+1))
371              q(i+1, jjp1, llm+1-l, 1) = 0.25*(sxs(i)+sxs(i+1))
372            END DO
373            q(1, 1, llm+1-l, 1) = q(iip1, 1, llm+1-l, 1)
374            q(1, jjp1, llm+1-l, 1) = q(iip1, jjp1, llm+1-l, 1)
375    
376          END IF
377    
378        END DO
379      END IF
380    
381      ! bouclage en longitude
382      DO iq = 0, 3
383        DO l = 1, llm
384          DO j = 1, jjp1
385            q(iip1, j, l, iq) = q(1, j, l, iq)
386          END DO
387        END DO
388      END DO
389    
390      DO l = 1, llm
391        DO j = 1, jjp1
392          DO i = 1, iip1
393            IF (q(i,j,l,0)<0.) THEN
394              q(i, j, l, 0) = 0.
395            END IF
396          END DO
397        END DO
398      END DO
399    
400      DO l = 1, llm
401        DO j = 1, jjp1
402          DO i = 1, iip1
403            IF (q(i,j,l,0)<qmin) PRINT *, 'apres pentes, s0(', i, ',', j, ',', l, &
404              ')=', q(i, j, l, 0)
405          END DO
406        END DO
407      END DO
408      RETURN
409    END SUBROUTINE pentes_ini
Legend:
Removed from v.3  
changed lines
  Added in v.82
  ViewVC Help
Powered by ViewVC 1.1.21