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	module dissip_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE dissip(vcov, ucov, teta, p, dv, du, dh)
8
9	! 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
14	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, cdivu, crot, &
18	cdivh
19	use gradiv2_m, only: gradiv2
20
21	! 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
26	! 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
32	INTEGER l, ij
33
34	!-----------------------------------------------------------------------
35
36	! Initializations:
37	te1dt = tetaudiv * dtdiss
38	te2dt = tetaurot * dtdiss
39	te3dt = tetah * dtdiss
40	du = 0.
41	dv = 0.
42	dh = 0.
43
44	! Calcul de la dissipation:
45
46	! Calcul de la partie grad (div) :
47
48	IF (lstardis) THEN
49	CALL gradiv2(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)
50	ELSE
51	CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu)
52	END IF
53
54	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
60	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
68	! calcul de la partie n X grad (rot) :
69
70	IF (lstardis) THEN
71	CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry, crot)
72	ELSE
73	CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry, crot)
74	END IF
75
76
77	DO l = 1, llm
78	DO ij = 1, iip1
79	grx(ij, l) = 0.
80	END DO
81
82	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
90	! calcul de la partie div (grad) :
91
92	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
99	CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh)
100	ELSE
101	CALL divgrad(llm, teta, niterh, gdx, cdivh)
102	END IF
103
104	forall (l = 1: llm) dh(:, l) = dh(:, l) - te3dt(l) * gdx(:, l)
105
106	END SUBROUTINE dissip
107
108	end module dissip_m