libf/dyn3d/dissip.f90

SUBROUTINE dissip(vcov, ucov, teta, p, dv, du, dh)

  ! 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

  IMPLICIT NONE

  ! Arguments: 
  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)

  ! Local: 
  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

  !-----------------------------------------------------------------------

  ! initialisations: 

  DO l = 1, llm
     te1dt(l) = tetaudiv(l)*dtdiss
     te2dt(l) = tetaurot(l)*dtdiss
     te3dt(l) = tetah(l)*dtdiss
  END DO
  du = 0.
  dv = 0.
  dh = 0.

  ! Calcul de la dissipation: 

  ! Calcul de la partie grad (div) : 

  IF (lstardis) THEN
     CALL gradiv2(llm, ucov, vcov, nitergdiv, gdx, gdy)
  ELSE
     CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy)
  END IF

  DO l = 1, llm
     DO ij = 1, iip1
        gdx(ij, l) = 0.
        gdx(ij+ip1jm, l) = 0.
     END DO

     DO ij = iip2, ip1jm
        du(ij, l) = du(ij, l) - te1dt(l)*gdx(ij, l)
     END DO
     DO ij = 1, ip1jm
        dv(ij, l) = dv(ij, l) - te1dt(l)*gdy(ij, l)
     END DO
  END DO

  ! calcul de la partie n X grad (rot) : 

  IF (lstardis) THEN
     CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry)
  ELSE
     CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry)
  END IF


  DO l = 1, llm
     DO ij = 1, iip1
        grx(ij, l) = 0.
     END DO

     DO ij = iip2, ip1jm
        du(ij, l) = du(ij, l) - te2dt(l)*grx(ij, l)
     END DO
     DO ij = 1, ip1jm
        dv(ij, l) = dv(ij, l) - te2dt(l)*gry(ij, l)
     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

END SUBROUTINE dissip
1	SUBROUTINE dissip(vcov, ucov, teta, p, dv, du, dh)
2
3	! From dyn3d/dissip.F, version 1.1.1.1 2004/05/19 12:53:05
4	! Avec nouveaux operateurs star : gradiv2 , divgrad2, nxgraro2
5	! Author: P. Le Van
6	! Objet : dissipation horizontale
7
8	USE dimens_m, ONLY : llm
9	USE paramet_m, ONLY : iip1, iip2, ip1jm, ip1jmp1, llmp1
10	USE comdissnew, ONLY : lstardis, nitergdiv, nitergrot, niterh
11	USE inidissip_m, ONLY : dtdiss, tetah, tetaudiv, tetaurot
12
13	IMPLICIT NONE
14
15	! Arguments:
16	REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm), teta(ip1jmp1, llm)
17	REAL, INTENT (IN) :: p(ip1jmp1, llmp1)
18	REAL dv(ip1jm, llm), du(ip1jmp1, llm), dh(ip1jmp1, llm)
19
20	! Local:
21	REAL gdx(ip1jmp1, llm), gdy(ip1jm, llm)
22	REAL grx(ip1jmp1, llm), gry(ip1jm, llm)
23	REAL te1dt(llm), te2dt(llm), te3dt(llm)
24	REAL deltapres(ip1jmp1, llm)
25
26	INTEGER l, ij
27
28	!-----------------------------------------------------------------------
29
30	! initialisations:
31
32	DO l = 1, llm
33	te1dt(l) = tetaudiv(l)*dtdiss
34	te2dt(l) = tetaurot(l)*dtdiss
35	te3dt(l) = tetah(l)*dtdiss
36	END DO
37	du = 0.
38	dv = 0.
39	dh = 0.
40
41	! Calcul de la dissipation:
42
43	! Calcul de la partie grad (div) :
44
45	IF (lstardis) THEN
46	CALL gradiv2(llm, ucov, vcov, nitergdiv, gdx, gdy)
47	ELSE
48	CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy)
49	END IF
50
51	DO l = 1, llm
52	DO ij = 1, iip1
53	gdx(ij, l) = 0.
54	gdx(ij+ip1jm, l) = 0.
55	END DO
56
57	DO ij = iip2, ip1jm
58	du(ij, l) = du(ij, l) - te1dt(l)*gdx(ij, l)
59	END DO
60	DO ij = 1, ip1jm
61	dv(ij, l) = dv(ij, l) - te1dt(l)*gdy(ij, l)
62	END DO
63	END DO
64
65	! calcul de la partie n X grad (rot) :
66
67	IF (lstardis) THEN
68	CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry)
69	ELSE
70	CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry)
71	END IF
72
73
74	DO l = 1, llm
75	DO ij = 1, iip1
76	grx(ij, l) = 0.
77	END DO
78
79	DO ij = iip2, ip1jm
80	du(ij, l) = du(ij, l) - te2dt(l)*grx(ij, l)
81	END DO
82	DO ij = 1, ip1jm
83	dv(ij, l) = dv(ij, l) - te2dt(l)*gry(ij, l)
84	END DO
85	END DO
86
87	! calcul de la partie div (grad) :
88
89	IF (lstardis) THEN
90	DO l = 1, llm
91	DO ij = 1, ip1jmp1
92	deltapres(ij, l) = amax1(0., p(ij, l)-p(ij, l+1))
93	END DO
94	END DO
95
96	CALL divgrad2(llm, teta, deltapres, niterh, gdx)
97	ELSE
98	CALL divgrad(llm, teta, niterh, gdx)
99	END IF
100
101	DO l = 1, llm
102	DO ij = 1, ip1jmp1
103	dh(ij, l) = dh(ij, l) - te3dt(l)*gdx(ij, l)
104	END DO
105	END DO
106
107	END SUBROUTINE dissip