libf/dyn3d/dissip.f90

module dissip_m

  IMPLICIT NONE

contains

  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 opérateurs star : gradiv2, divgrad2, nxgraro2
    ! Author: P. Le Van
    ! Objet : dissipation horizontale

    USE dimens_m, ONLY : iim, jjm, llm
    USE paramet_m, ONLY : iip1, iip2, ip1jmp1, llmp1
    USE comdissnew, ONLY : lstardis, nitergdiv, nitergrot, niterh
    USE inidissip_m, ONLY : dtdiss, tetah, tetaudiv, tetaurot

    ! Arguments:
    REAL vcov((iim + 1) * jjm, llm), ucov(ip1jmp1, llm), teta(ip1jmp1, llm)
    REAL, INTENT (IN) :: p(ip1jmp1, llmp1)
    REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm), dh(ip1jmp1, llm)

    ! Local:
    REAL gdx(ip1jmp1, llm), gdy((iim + 1) * jjm, llm)
    REAL grx(ip1jmp1, llm), gry((iim + 1) * jjm, llm)
    REAL te1dt(llm), te2dt(llm), te3dt(llm)
    REAL deltapres(ip1jmp1, llm)

    INTEGER l, ij

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

    ! Initializations:
    te1dt = tetaudiv * dtdiss
    te2dt = tetaurot * dtdiss
    te3dt = tetah * dtdiss
    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+(iim + 1) * jjm, l) = 0.
       END DO

       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

    ! 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, (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

    ! calcul de la partie div (grad) :

    IF (lstardis) THEN
       DO l = 1, llm
          DO ij = 1, ip1jmp1
             deltapres(ij, l) = max(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

    forall (l = 1: llm) dh(:, l) = dh(:, l) - te3dt(l) * gdx(:, l)

  END SUBROUTINE dissip

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