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