/[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 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/dyn3d/Dissipation/dissip.f revision 279 by guez, Fri Jul 20 14:30:23 2018 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  
       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   -------------  
   
       include "comdissipn.h"  
   
 c   Arguments:  
 c   ----------  
   
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)  
       REAL  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  
   
 c   calcul de la partie   div ( grad ):  
 c   -----------------------------------  
   
           
       IF(lstardis) THEN  
   
        DO l = 1, llm  
           DO ij = 1, ip1jmp1  
             deltapres(ij,l) = AMAX1( 0.,  p(ij,l) - p(ij,l+1) )  
           ENDDO  
        ENDDO  
   
          CALL divgrad2( llm,teta, deltapres  ,niterh, gdx )  
       ELSE  
          CALL divgrad ( llm,teta, niterh, gdx        )  
       ENDIF  
   
       DO l = 1,llm  
          DO ij = 1,ip1jmp1  
             dh( ij,l ) = dh( ij,l ) - te3dt(l) * gdx( ij,l )  
          ENDDO  
       ENDDO  
2    
3        RETURN    IMPLICIT NONE
4        END  
5    contains
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        ! Author: P. Le Van
11        
12        ! Objet : calcul de la dissipation horizontale. Avec op\'erateurs
13        ! star : gradiv2, divgrad2, nxgraro2.
14    
15        use nr_util, only: assert
16    
17        USE comdissnew, ONLY: nitergdiv, nitergrot, niterh
18        USE dimensions, ONLY: iim, jjm, llm
19        use divgrad2_m, only: divgrad2
20        use gradiv2_m, only: gradiv2
21        USE inidissip_m, ONLY: dtdiss, tetah, tetaudiv, tetaurot, cdivu, crot, cdivh
22        use nxgraro2_m, only: nxgraro2
23    
24        REAL, intent(in):: vcov(:, :, :) ! (iim + 1, jjm, llm)
25        REAL, intent(in):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm)
26        REAL, intent(in):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
27        REAL, INTENT(IN):: p(:, :, :) ! (iim + 1, jjm + 1, llm + 1)
28        REAL, intent(out):: dv(:, :, :) ! (iim + 1, jjm, llm)
29        REAL, intent(out):: du(:, :, :) ! (iim + 1, jjm + 1, llm)
30        REAL, intent(out):: dh(:, :, :) ! (iim + 1, jjm + 1, llm)
31    
32        ! Local:
33        REAL gdx(iim + 1, jjm + 1, llm), gdy(iim + 1, jjm, llm)
34        REAL tedt(llm)
35        REAL deltapres(iim + 1, jjm + 1, llm)
36        INTEGER l
37    
38        !-----------------------------------------------------------------------
39    
40        call assert((/size(vcov, 1), size(ucov, 1), size(teta, 1), size(p, 1), &
41             size(dv, 1), size(du, 1), size(dh, 1)/) == iim + 1, "dissip iim")
42        call assert((/size(vcov, 2), size(ucov, 2) - 1, size(teta, 2) - 1, &
43             size(p, 2) - 1, size(dv, 2), size(du, 2) - 1, size(dh, 2) - 1/) &
44             == jjm, "dissip jjm")
45        call assert((/size(vcov, 3), size(ucov, 3), size(teta, 3), size(p, 3) - 1, &
46             size(dv, 3), size(du, 3), size(dh, 3)/) == llm, "dissip llm")
47    
48        du(:, 1, :) = 0.
49        du(:, jjm + 1, :) = 0.
50    
51        ! Calcul de la partie grad(div) :
52        CALL gradiv2(ucov, vcov, nitergdiv, gdx, gdy, cdivu)
53        tedt = tetaudiv * dtdiss
54        forall (l = 1: llm)
55           du(:, 2: jjm, l) = - tedt(l) * gdx(:, 2: jjm, l)
56           dv(:, :, l) = - tedt(l) * gdy(:, :, l)
57        END forall
58    
59        ! Calcul de la partie n \wedge grad(rot) :
60        CALL nxgraro2(ucov, vcov, nitergrot, gdx, gdy, crot)
61        tedt = tetaurot * dtdiss
62        forall (l = 1: llm)
63           du(:, 2: jjm, l) = du(:, 2: jjm, l) - tedt(l) * gdx(:, 2: jjm, l)
64           dv(:, :, l) = dv(:, :, l) - tedt(l) * gdy(:, :, l)
65        END forall
66    
67        ! calcul de la partie div(grad) :
68        forall (l = 1: llm) &
69             deltapres(:, :, l) = max(0., p(:, :, l) - p(:, :, l + 1))
70        CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh)
71        forall (l = 1: llm) dh(:, :, l) = - tetah(l) * dtdiss * gdx(:, :, l)
72    
73      END SUBROUTINE dissip
74    
75    end module dissip_m
Legend:
Removed from v.3  
changed lines
  Added in v.279
  ViewVC Help
Powered by ViewVC 1.1.21