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, cdivu, crot, &
         cdivh
    use gradiv2_m, only: gradiv2

    ! 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, cdivu)
    ELSE
       CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)
    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, crot)
    ELSE
       CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry, crot)
    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, cdivh)
    ELSE
       CALL divgrad(llm, teta, niterh, gdx, cdivh)
    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	guez	54	USE inidissip_m, ONLY : dtdiss, tetah, tetaudiv, tetaurot, cdivu, crot, &
18			cdivh
19			use gradiv2_m, only: gradiv2
20	guez	3
21	guez	47	! Arguments:
22			REAL vcov((iim + 1) * jjm, llm), ucov(ip1jmp1, llm), teta(ip1jmp1, llm)
23			REAL, INTENT (IN) :: p(ip1jmp1, llmp1)
24			REAL dv((iim + 1) * jjm, llm), du(ip1jmp1, llm), dh(ip1jmp1, llm)
25	guez	3
26	guez	47	! Local:
27			REAL gdx(ip1jmp1, llm), gdy((iim + 1) * jjm, llm)
28			REAL grx(ip1jmp1, llm), gry((iim + 1) * jjm, llm)
29			REAL te1dt(llm), te2dt(llm), te3dt(llm)
30			REAL deltapres(ip1jmp1, llm)
31	guez	3
32	guez	47	INTEGER l, ij
33	guez	3
34	guez	47	!-----------------------------------------------------------------------
35	guez	3
36	guez	47	! Initializations:
37			te1dt = tetaudiv * dtdiss
38			te2dt = tetaurot * dtdiss
39			te3dt = tetah * dtdiss
40			du = 0.
41			dv = 0.
42			dh = 0.
43	guez	3
44	guez	47	! Calcul de la dissipation:
45	guez	3
46	guez	47	! Calcul de la partie grad (div) :
47	guez	3
48	guez	47	IF (lstardis) THEN
49	guez	54	CALL gradiv2(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)
50	guez	47	ELSE
51	guez	54	CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)
52	guez	47	END IF
53	guez	3
54	guez	47	DO l = 1, llm
55			DO ij = 1, iip1
56			gdx(ij, l) = 0.
57			gdx(ij+(iim + 1) * jjm, l) = 0.
58			END DO
59	guez	3
60	guez	47	DO ij = iip2, (iim + 1) * jjm
61			du(ij, l) = du(ij, l) - te1dt(l) * gdx(ij, l)
62			END DO
63			DO ij = 1, (iim + 1) * jjm
64			dv(ij, l) = dv(ij, l) - te1dt(l) * gdy(ij, l)
65			END DO
66			END DO
67	guez	3
68	guez	47	! calcul de la partie n X grad (rot) :
69	guez	3
70	guez	47	IF (lstardis) THEN
71	guez	54	CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry, crot)
72	guez	47	ELSE
73	guez	54	CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry, crot)
74	guez	47	END IF
75	guez	3
76
77	guez	47	DO l = 1, llm
78			DO ij = 1, iip1
79			grx(ij, l) = 0.
80			END DO
81	guez	3
82	guez	47	DO ij = iip2, (iim + 1) * jjm
83			du(ij, l) = du(ij, l) - te2dt(l) * grx(ij, l)
84			END DO
85			DO ij = 1, (iim + 1) * jjm
86			dv(ij, l) = dv(ij, l) - te2dt(l) * gry(ij, l)
87			END DO
88			END DO
89	guez	3
90	guez	47	! calcul de la partie div (grad) :
91	guez	3
92	guez	47	IF (lstardis) THEN
93			DO l = 1, llm
94			DO ij = 1, ip1jmp1
95			deltapres(ij, l) = max(0., p(ij, l) - p(ij, l + 1))
96			END DO
97			END DO
98	guez	3
99	guez	54	CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh)
100	guez	47	ELSE
101	guez	54	CALL divgrad(llm, teta, niterh, gdx, cdivh)
102	guez	47	END IF
103	guez	3
104	guez	47	forall (l = 1: llm) dh(:, l) = dh(:, l) - te3dt(l) * gdx(:, l)
105
106			END SUBROUTINE dissip
107
108			end module dissip_m