/[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.f90 revision 54 by guez, Tue Dec 6 15:07:04 2011 UTC trunk/libf/dyn3d/Dissipation/dissip.f90 revision 57 by guez, Mon Jan 30 12:54:02 2012 UTC
# Line 7  contains Line 7  contains
7    SUBROUTINE dissip(vcov, ucov, teta, p, dv, du, dh)    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      ! From dyn3d/dissip.F, version 1.1.1.1 2004/05/19 12:53:05
     ! Avec nouveaux opérateurs star : gradiv2, divgrad2, nxgraro2  
10      ! Author: P. Le Van      ! Author: P. Le Van
11      ! Objet : dissipation horizontale      ! Objet : calcul de la dissipation horizontale
12        ! Avec opérateurs star : gradiv2, divgrad2, nxgraro2
13    
14      USE dimens_m, ONLY : iim, jjm, llm      USE dimens_m, ONLY: iim, jjm, llm
15      USE paramet_m, ONLY : iip1, iip2, ip1jmp1, llmp1      USE comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh
16      USE comdissnew, ONLY : lstardis, nitergdiv, nitergrot, niterh      USE inidissip_m, ONLY: dtdiss, tetah, tetaudiv, tetaurot, cdivu, crot, cdivh
     USE inidissip_m, ONLY : dtdiss, tetah, tetaudiv, tetaurot, cdivu, crot, &  
          cdivh  
17      use gradiv2_m, only: gradiv2      use gradiv2_m, only: gradiv2
18        use nr_util, only: assert
19    
20      ! Arguments:      REAL, intent(in):: vcov(:, :, :) ! (iim + 1, jjm, llm)
21      REAL vcov((iim + 1) * jjm, llm), ucov(ip1jmp1, llm), teta(ip1jmp1, llm)      REAL, intent(in):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm)
22      REAL, INTENT (IN) :: p(ip1jmp1, llmp1)      REAL, intent(in):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
23      REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm), dh(ip1jmp1, llm)      REAL, INTENT(IN):: p(:, :, :) ! (iim + 1, jjm + 1, llm + 1)
24        REAL, intent(out):: dv(:, :, :) ! (iim + 1, jjm, llm)
25        REAL, intent(out):: du(:, :, :) ! (iim + 1, jjm + 1, llm)
26        REAL, intent(out):: dh(:, :, :) ! (iim + 1, jjm + 1, llm)
27    
28      ! Local:      ! Local:
29      REAL gdx(ip1jmp1, llm), gdy((iim + 1) * jjm, llm)      REAL gdx(iim + 1, jjm + 1, llm), gdy(iim + 1, jjm, llm)
30      REAL grx(ip1jmp1, llm), gry((iim + 1) * jjm, llm)      REAL grx(iim + 1, jjm + 1, llm), gry(iim + 1, jjm, llm)
31      REAL te1dt(llm), te2dt(llm), te3dt(llm)      REAL tedt(llm)
32      REAL deltapres(ip1jmp1, llm)      REAL deltapres(iim + 1, jjm + 1, llm)
33        INTEGER l
     INTEGER l, ij  
34    
35      !-----------------------------------------------------------------------      !-----------------------------------------------------------------------
36    
37      ! Initializations:      call assert((/size(vcov, 1), size(ucov, 1), size(teta, 1), size(p, 1), &
38      te1dt = tetaudiv * dtdiss           size(dv, 1), size(du, 1), size(dh, 1)/) == iim + 1, "dissip iim")
39      te2dt = tetaurot * dtdiss      call assert((/size(vcov, 2), size(ucov, 2) - 1, size(teta, 2) - 1, &
40      te3dt = tetah * dtdiss           size(p, 2) - 1, size(dv, 2), size(du, 2) - 1, size(dh, 2) - 1/) &
41      du = 0.           == jjm, "dissip jjm")
42      dv = 0.      call assert((/size(vcov, 3), size(ucov, 3), size(teta, 3), size(p, 3) - 1, &
43      dh = 0.           size(dv, 3), size(du, 3), size(dh, 3)/) == llm, "dissip llm")
44    
45      ! Calcul de la dissipation:      du(:, 1, :) = 0.
46        du(:, jjm + 1, :) = 0.
47    
48      ! Calcul de la partie grad (div) :      ! Calcul de la partie grad (div) :
49    
50      IF (lstardis) THEN      IF (lstardis) THEN
51         CALL gradiv2(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)         CALL gradiv2(ucov, vcov, nitergdiv, gdx, gdy, cdivu)
52      ELSE      ELSE
53         CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)         CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)
54      END IF      END IF
55    
56      DO l = 1, llm      tedt = tetaudiv * dtdiss
57         DO ij = 1, iip1      forall (l = 1: llm)
58            gdx(ij, l) = 0.         du(:, 2: jjm, l) = - tedt(l) * gdx(:, 2: jjm, l)
59            gdx(ij+(iim + 1) * jjm, l) = 0.         dv(:, :, l) = - tedt(l) * gdy(:, :, l)
60         END DO      END forall
   
        DO ij = iip2, (iim + 1) * jjm  
           du(ij, l) = du(ij, l) - te1dt(l) * gdx(ij, l)  
        END DO  
        DO ij = 1, (iim + 1) * jjm  
           dv(ij, l) = dv(ij, l) - te1dt(l) * gdy(ij, l)  
        END DO  
     END DO  
61    
62      ! calcul de la partie n X grad (rot) :      ! Calcul de la partie n X grad (rot) :
63    
64      IF (lstardis) THEN      IF (lstardis) THEN
65         CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry, crot)         CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry, crot)
# Line 73  contains Line 67  contains
67         CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry, crot)         CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry, crot)
68      END IF      END IF
69    
70        tedt = tetaurot * dtdiss
71      DO l = 1, llm      forall (l = 1: llm)
72         DO ij = 1, iip1         du(:, 2: jjm, l) = du(:, 2: jjm, l) - tedt(l) * grx(:, 2: jjm, l)
73            grx(ij, l) = 0.         dv(:, :, l) = dv(:, :, l) - tedt(l) * gry(:, :, l)
74         END DO      END forall
   
        DO ij = iip2, (iim + 1) * jjm  
           du(ij, l) = du(ij, l) - te2dt(l) * grx(ij, l)  
        END DO  
        DO ij = 1, (iim + 1) * jjm  
           dv(ij, l) = dv(ij, l) - te2dt(l) * gry(ij, l)  
        END DO  
     END DO  
75    
76      ! calcul de la partie div (grad) :      ! calcul de la partie div (grad) :
77    
78      IF (lstardis) THEN      IF (lstardis) THEN
79         DO l = 1, llm         forall (l = 1: llm) &
80            DO ij = 1, ip1jmp1              deltapres(:, :, l) = max(0., p(:, :, l) - p(:, :, l + 1))
              deltapres(ij, l) = max(0., p(ij, l) - p(ij, l + 1))  
           END DO  
        END DO  
   
81         CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh)         CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh)
82      ELSE      ELSE
83         CALL divgrad(llm, teta, niterh, gdx, cdivh)         CALL divgrad(llm, teta, niterh, gdx, cdivh)
84      END IF      END IF
85    
86      forall (l = 1: llm) dh(:, l) = dh(:, l) - te3dt(l) * gdx(:, l)      forall (l = 1: llm) dh(:, :, l) = - tetah(l) * dtdiss * gdx(:, :, l)
87    
88    END SUBROUTINE dissip    END SUBROUTINE dissip
89    
Legend:
Removed from v.54  
changed lines
  Added in v.57
  ViewVC Help
Powered by ViewVC 1.1.21