dyn3d/Dissipation/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
    ! Author: P. Le Van
    
    ! Objet : calcul de la dissipation horizontale. Avec op\'erateurs
    ! star : gradiv2, divgrad2, nxgraro2.

    use nr_util, only: assert

    USE comdissnew, ONLY: nitergdiv, nitergrot, niterh
    USE dimensions, ONLY: iim, jjm, llm
    use divgrad2_m, only: divgrad2
    use gradiv2_m, only: gradiv2
    USE inidissip_m, ONLY: dtdiss, tetah, tetaudiv, tetaurot, cdivu, crot, cdivh
    use nxgraro2_m, only: nxgraro2

    REAL, intent(in):: vcov(:, :, :) ! (iim + 1, jjm, llm)
    REAL, intent(in):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm)
    REAL, intent(in):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
    REAL, INTENT(IN):: p(:, :, :) ! (iim + 1, jjm + 1, llm + 1)
    REAL, intent(out):: dv(:, :, :) ! (iim + 1, jjm, llm)
    REAL, intent(out):: du(:, :, :) ! (iim + 1, jjm + 1, llm)
    REAL, intent(out):: dh(:, :, :) ! (iim + 1, jjm + 1, llm)

    ! Local:
    REAL gdx(iim + 1, jjm + 1, llm), gdy(iim + 1, jjm, llm)
    REAL tedt(llm)
    REAL deltapres(iim + 1, jjm + 1, llm)
    INTEGER l

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

    call assert((/size(vcov, 1), size(ucov, 1), size(teta, 1), size(p, 1), &
         size(dv, 1), size(du, 1), size(dh, 1)/) == iim + 1, "dissip iim")
    call assert((/size(vcov, 2), size(ucov, 2) - 1, size(teta, 2) - 1, &
         size(p, 2) - 1, size(dv, 2), size(du, 2) - 1, size(dh, 2) - 1/) &
         == jjm, "dissip jjm")
    call assert((/size(vcov, 3), size(ucov, 3), size(teta, 3), size(p, 3) - 1, &
         size(dv, 3), size(du, 3), size(dh, 3)/) == llm, "dissip llm")

    du(:, 1, :) = 0.
    du(:, jjm + 1, :) = 0.

    ! Calcul de la partie grad(div) :
    CALL gradiv2(ucov, vcov, nitergdiv, gdx, gdy, cdivu)
    tedt = tetaudiv * dtdiss
    forall (l = 1: llm)
       du(:, 2: jjm, l) = - tedt(l) * gdx(:, 2: jjm, l)
       dv(:, :, l) = - tedt(l) * gdy(:, :, l)
    END forall

    ! Calcul de la partie n \wedge grad(rot) :
    CALL nxgraro2(ucov, vcov, nitergrot, gdx, gdy, crot)
    tedt = tetaurot * dtdiss
    forall (l = 1: llm)
       du(:, 2: jjm, l) = du(:, 2: jjm, l) - tedt(l) * gdx(:, 2: jjm, l)
       dv(:, :, l) = dv(:, :, l) - tedt(l) * gdy(:, :, l)
    END forall

    ! calcul de la partie div(grad) :
    forall (l = 1: llm) &
         deltapres(:, :, l) = max(0., p(:, :, l) - p(:, :, l + 1))
    CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh)
    forall (l = 1: llm) dh(:, :, l) = - tetah(l) * dtdiss * 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			! Author: P. Le Van
11	guez	279
12			! Objet : calcul de la dissipation horizontale. Avec op\'erateurs
13			! star : gradiv2, divgrad2, nxgraro2.
14	guez	3
15	guez	279	use nr_util, only: assert
16
17	guez	65	USE comdissnew, ONLY: nitergdiv, nitergrot, niterh
18	guez	265	USE dimensions, ONLY: iim, jjm, llm
19	guez	65	use divgrad2_m, only: divgrad2
20			use gradiv2_m, only: gradiv2
21	guez	55	USE inidissip_m, ONLY: dtdiss, tetah, tetaudiv, tetaurot, cdivu, crot, cdivh
22	guez	65	use nxgraro2_m, only: nxgraro2
23	guez	3
24	guez	55	REAL, intent(in):: vcov(:, :, :) ! (iim + 1, jjm, llm)
25			REAL, intent(in):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm)
26	guez	56	REAL, intent(in):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
27			REAL, INTENT(IN):: p(:, :, :) ! (iim + 1, jjm + 1, llm + 1)
28	guez	55	REAL, intent(out):: dv(:, :, :) ! (iim + 1, jjm, llm)
29			REAL, intent(out):: du(:, :, :) ! (iim + 1, jjm + 1, llm)
30			REAL, intent(out):: dh(:, :, :) ! (iim + 1, jjm + 1, llm)
31	guez	3
32	guez	47	! Local:
33	guez	55	REAL gdx(iim + 1, jjm + 1, llm), gdy(iim + 1, jjm, llm)
34	guez	56	REAL tedt(llm)
35			REAL deltapres(iim + 1, jjm + 1, llm)
36	guez	55	INTEGER l
37	guez	3
38	guez	47	!-----------------------------------------------------------------------
39	guez	3
40	guez	56	call assert((/size(vcov, 1), size(ucov, 1), size(teta, 1), size(p, 1), &
41			size(dv, 1), size(du, 1), size(dh, 1)/) == iim + 1, "dissip iim")
42			call assert((/size(vcov, 2), size(ucov, 2) - 1, size(teta, 2) - 1, &
43			size(p, 2) - 1, size(dv, 2), size(du, 2) - 1, size(dh, 2) - 1/) &
44			== jjm, "dissip jjm")
45			call assert((/size(vcov, 3), size(ucov, 3), size(teta, 3), size(p, 3) - 1, &
46			size(dv, 3), size(du, 3), size(dh, 3)/) == llm, "dissip llm")
47	guez	3
48	guez	56	du(:, 1, :) = 0.
49			du(:, jjm + 1, :) = 0.
50	guez	3
51	guez	65	! Calcul de la partie grad(div) :
52	guez	64	CALL gradiv2(ucov, vcov, nitergdiv, gdx, gdy, cdivu)
53	guez	56	tedt = tetaudiv * dtdiss
54	guez	55	forall (l = 1: llm)
55	guez	56	du(:, 2: jjm, l) = - tedt(l) * gdx(:, 2: jjm, l)
56			dv(:, :, l) = - tedt(l) * gdy(:, :, l)
57	guez	55	END forall
58	guez	3
59	guez	279	! Calcul de la partie n \wedge grad(rot) :
60	guez	65	CALL nxgraro2(ucov, vcov, nitergrot, gdx, gdy, crot)
61	guez	56	tedt = tetaurot * dtdiss
62	guez	55	forall (l = 1: llm)
63	guez	65	du(:, 2: jjm, l) = du(:, 2: jjm, l) - tedt(l) * gdx(:, 2: jjm, l)
64			dv(:, :, l) = dv(:, :, l) - tedt(l) * gdy(:, :, l)
65	guez	55	END forall
66	guez	3
67	guez	65	! calcul de la partie div(grad) :
68	guez	64	forall (l = 1: llm) &
69			deltapres(:, :, l) = max(0., p(:, :, l) - p(:, :, l + 1))
70			CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh)
71	guez	56	forall (l = 1: llm) dh(:, :, l) = - tetah(l) * dtdiss * gdx(:, :, l)
72	guez	47
73			END SUBROUTINE dissip
74
75			end module dissip_m