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	guez	44	SUBROUTINE dissip(vcov, ucov, teta, p, dv, du, dh)
2	guez	3
3	guez	44	! 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	guez	3
8	guez	26	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	guez	3
13	guez	26	IMPLICIT NONE
14	guez	3
15	guez	44	! 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	guez	3
20	guez	44	! 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	guez	3
26	guez	44	INTEGER l, ij
27	guez	3
28	guez	26	!-----------------------------------------------------------------------
29	guez	3
30	guez	44	! initialisations:
31	guez	3
32	guez	26	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	guez	3
41	guez	44	! Calcul de la dissipation:
42	guez	3
43	guez	44	! Calcul de la partie grad (div) :
44	guez	3
45	guez	26	IF (lstardis) THEN
46	guez	44	CALL gradiv2(llm, ucov, vcov, nitergdiv, gdx, gdy)
47	guez	26	ELSE
48	guez	44	CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy)
49	guez	26	END IF
50	guez	3
51	guez	26	DO l = 1, llm
52			DO ij = 1, iip1
53	guez	44	gdx(ij, l) = 0.
54			gdx(ij+ip1jm, l) = 0.
55	guez	26	END DO
56	guez	3
57	guez	26	DO ij = iip2, ip1jm
58	guez	44	du(ij, l) = du(ij, l) - te1dt(l)*gdx(ij, l)
59	guez	26	END DO
60			DO ij = 1, ip1jm
61	guez	44	dv(ij, l) = dv(ij, l) - te1dt(l)*gdy(ij, l)
62	guez	26	END DO
63			END DO
64	guez	3
65	guez	44	! calcul de la partie n X grad (rot) :
66	guez	3
67	guez	26	IF (lstardis) THEN
68	guez	44	CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry)
69	guez	26	ELSE
70	guez	44	CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry)
71	guez	26	END IF
72	guez	3
73
74	guez	26	DO l = 1, llm
75			DO ij = 1, iip1
76	guez	44	grx(ij, l) = 0.
77	guez	26	END DO
78	guez	3
79	guez	26	DO ij = iip2, ip1jm
80	guez	44	du(ij, l) = du(ij, l) - te2dt(l)*grx(ij, l)
81	guez	26	END DO
82			DO ij = 1, ip1jm
83	guez	44	dv(ij, l) = dv(ij, l) - te2dt(l)*gry(ij, l)
84	guez	26	END DO
85			END DO
86	guez	3
87	guez	44	! calcul de la partie div (grad) :
88	guez	3
89	guez	26	IF (lstardis) THEN
90			DO l = 1, llm
91			DO ij = 1, ip1jmp1
92	guez	44	deltapres(ij, l) = amax1(0., p(ij, l)-p(ij, l+1))
93	guez	26	END DO
94			END DO
95	guez	3
96	guez	44	CALL divgrad2(llm, teta, deltapres, niterh, gdx)
97	guez	26	ELSE
98	guez	44	CALL divgrad(llm, teta, niterh, gdx)
99	guez	26	END IF
100	guez	3
101	guez	26	DO l = 1, llm
102			DO ij = 1, ip1jmp1
103	guez	44	dh(ij, l) = dh(ij, l) - te3dt(l)*gdx(ij, l)
104	guez	26	END DO
105			END DO
106	guez	3
107	guez	26	END SUBROUTINE dissip