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