/[lmdze]/trunk/dyn3d/Dissipation/dissip.f
ViewVC logotype

Diff of /trunk/dyn3d/Dissipation/dissip.f

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

trunk/libf/dyn3d/dissip.f revision 25 by guez, Fri Mar 5 16:43:45 2010 UTC trunk/libf/dyn3d/dissip.f90 revision 54 by guez, Tue Dec 6 15:07:04 2011 UTC
# Line 1  Line 1 
1  !  module dissip_m
 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/dissip.F,v 1.1.1.1 2004/05/19 12:53:05 lmdzadmin Exp $  
 !  
       SUBROUTINE dissip( vcov,ucov,teta,p, dv,du,dh )  
 c  
       use dimens_m  
       use paramet_m  
       use comconst  
       use comdissnew  
       use comgeom  
             use comdissipn  
       IMPLICIT NONE  
   
   
 c ..  Avec nouveaux operateurs star :  gradiv2 , divgrad2, nxgraro2  ...  
 c                                 (  10/01/98  )  
   
 c=======================================================================  
 c  
 c   Auteur:  P. Le Van  
 c   -------  
 c  
 c   Objet:  
 c   ------  
 c  
 c   Dissipation horizontale  
 c  
 c=======================================================================  
 c-----------------------------------------------------------------------  
 c   Declarations:  
 c   -------------  
   
   
 c   Arguments:  
 c   ----------  
   
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)  
       REAL, intent(in):: p( ip1jmp1,llmp1 )  
       REAL dv(ip1jm,llm),du(ip1jmp1,llm),dh(ip1jmp1,llm)  
   
 c   Local:  
 c   ------  
   
       REAL gdx(ip1jmp1,llm),gdy(ip1jm,llm)  
       REAL grx(ip1jmp1,llm),gry(ip1jm,llm)  
       REAL te1dt(llm),te2dt(llm),te3dt(llm)  
       REAL deltapres(ip1jmp1,llm)  
   
       INTEGER l,ij  
   
       REAL  SSUM  
   
 c-----------------------------------------------------------------------  
 c   initialisations:  
 c   ----------------  
   
       DO l=1,llm  
          te1dt(l) = tetaudiv(l) * dtdiss  
          te2dt(l) = tetaurot(l) * dtdiss  
          te3dt(l) = tetah(l)    * dtdiss  
       ENDDO  
       du=0.  
       dv=0.  
       dh=0.  
   
 c-----------------------------------------------------------------------  
 c   Calcul de la dissipation:  
 c   -------------------------  
   
 c   Calcul de la partie   grad  ( div ) :  
 c   -------------------------------------  
   
   
       IF(lstardis) THEN  
          CALL gradiv2( llm,ucov,vcov,nitergdiv,gdx,gdy )  
       ELSE  
          CALL gradiv ( llm,ucov,vcov,nitergdiv,gdx,gdy )  
       ENDIF  
   
       DO l=1,llm  
   
          DO ij = 1, iip1  
             gdx(     ij ,l) = 0.  
             gdx(ij+ip1jm,l) = 0.  
          ENDDO  
   
          DO ij = iip2,ip1jm  
             du(ij,l) = du(ij,l) - te1dt(l) *gdx(ij,l)  
          ENDDO  
          DO ij = 1,ip1jm  
             dv(ij,l) = dv(ij,l) - te1dt(l) *gdy(ij,l)  
          ENDDO  
   
        ENDDO  
   
 c   calcul de la partie   n X grad ( rot ):  
 c   ---------------------------------------  
   
       IF(lstardis) THEN  
          CALL nxgraro2( llm,ucov, vcov, nitergrot,grx,gry )  
       ELSE  
          CALL nxgrarot( llm,ucov, vcov, nitergrot,grx,gry )  
       ENDIF  
   
   
       DO l=1,llm  
          DO ij = 1, iip1  
             grx(ij,l) = 0.  
          ENDDO  
   
          DO ij = iip2,ip1jm  
             du(ij,l) = du(ij,l) - te2dt(l) * grx(ij,l)  
          ENDDO  
          DO ij =  1, ip1jm  
             dv(ij,l) = dv(ij,l) - te2dt(l) * gry(ij,l)  
          ENDDO  
       ENDDO  
2    
3  c   calcul de la partie   div ( grad ):    IMPLICIT NONE
 c   -----------------------------------  
4    
5            contains
       IF(lstardis) THEN  
6    
7      SUBROUTINE dissip(vcov, ucov, teta, p, dv, du, dh)
8    
9        ! From dyn3d/dissip.F, version 1.1.1.1 2004/05/19 12:53:05
10        ! Avec nouveaux opérateurs star : gradiv2, divgrad2, nxgraro2
11        ! Author: P. Le Van
12        ! Objet : dissipation horizontale
13    
14        USE dimens_m, ONLY : iim, jjm, llm
15        USE paramet_m, ONLY : iip1, iip2, ip1jmp1, llmp1
16        USE comdissnew, ONLY : lstardis, nitergdiv, nitergrot, niterh
17        USE inidissip_m, ONLY : dtdiss, tetah, tetaudiv, tetaurot, cdivu, crot, &
18             cdivh
19        use gradiv2_m, only: gradiv2
20    
21        ! Arguments:
22        REAL vcov((iim + 1) * jjm, llm), ucov(ip1jmp1, llm), teta(ip1jmp1, llm)
23        REAL, INTENT (IN) :: p(ip1jmp1, llmp1)
24        REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm), dh(ip1jmp1, llm)
25    
26        ! Local:
27        REAL gdx(ip1jmp1, llm), gdy((iim + 1) * jjm, llm)
28        REAL grx(ip1jmp1, llm), gry((iim + 1) * jjm, llm)
29        REAL te1dt(llm), te2dt(llm), te3dt(llm)
30        REAL deltapres(ip1jmp1, llm)
31    
32        INTEGER l, ij
33    
34        !-----------------------------------------------------------------------
35    
36        ! Initializations:
37        te1dt = tetaudiv * dtdiss
38        te2dt = tetaurot * dtdiss
39        te3dt = tetah * dtdiss
40        du = 0.
41        dv = 0.
42        dh = 0.
43    
44        ! Calcul de la dissipation:
45    
46        ! Calcul de la partie grad (div) :
47    
48        IF (lstardis) THEN
49           CALL gradiv2(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)
50        ELSE
51           CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)
52        END IF
53    
54        DO l = 1, llm
55           DO ij = 1, iip1
56              gdx(ij, l) = 0.
57              gdx(ij+(iim + 1) * jjm, l) = 0.
58           END DO
59    
60           DO ij = iip2, (iim + 1) * jjm
61              du(ij, l) = du(ij, l) - te1dt(l) * gdx(ij, l)
62           END DO
63           DO ij = 1, (iim + 1) * jjm
64              dv(ij, l) = dv(ij, l) - te1dt(l) * gdy(ij, l)
65           END DO
66        END DO
67    
68        ! calcul de la partie n X grad (rot) :
69    
70        IF (lstardis) THEN
71           CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry, crot)
72        ELSE
73           CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry, crot)
74        END IF
75    
76    
77        DO l = 1, llm
78           DO ij = 1, iip1
79              grx(ij, l) = 0.
80           END DO
81    
82           DO ij = iip2, (iim + 1) * jjm
83              du(ij, l) = du(ij, l) - te2dt(l) * grx(ij, l)
84           END DO
85           DO ij = 1, (iim + 1) * jjm
86              dv(ij, l) = dv(ij, l) - te2dt(l) * gry(ij, l)
87           END DO
88        END DO
89    
90        ! calcul de la partie div (grad) :
91    
92        IF (lstardis) THEN
93         DO l = 1, llm         DO l = 1, llm
94            DO ij = 1, ip1jmp1            DO ij = 1, ip1jmp1
95              deltapres(ij,l) = AMAX1( 0.,  p(ij,l) - p(ij,l+1) )               deltapres(ij, l) = max(0., p(ij, l) - p(ij, l + 1))
96            ENDDO            END DO
97         ENDDO         END DO
98    
99           CALL divgrad2( llm,teta, deltapres  ,niterh, gdx )         CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh)
100        ELSE      ELSE
101           CALL divgrad ( llm,teta, niterh, gdx        )         CALL divgrad(llm, teta, niterh, gdx, cdivh)
102        ENDIF      END IF
103    
104        DO l = 1,llm      forall (l = 1: llm) dh(:, l) = dh(:, l) - te3dt(l) * gdx(:, l)
105           DO ij = 1,ip1jmp1  
106              dh( ij,l ) = dh( ij,l ) - te3dt(l) * gdx( ij,l )    END SUBROUTINE dissip
          ENDDO  
       ENDDO  
107    
108        RETURN  end module dissip_m
       END  
Legend:
Removed from v.25  
changed lines
  Added in v.54
  ViewVC Help
Powered by ViewVC 1.1.21