trunk/dyn3d/advect.f

module advect_m

  IMPLICIT NONE

contains

  SUBROUTINE advect(ucov, vcov, teta, w, massebx, masseby, du, dv, dteta, &
       conser)

    ! From dyn3d/advect.F, version 1.1.1.1 2004/05/19 12:53:06
    ! Authors: P. Le Van , F. Hourdin
    ! Objet : calcul des termes d'advection verticale pour u, v, teta.
    ! Ces termes sont ajoutés à du, dv, dteta.

    USE dimens_m, ONLY : iim, llm
    USE paramet_m, ONLY : iip1, iip2, ip1jm, ip1jmp1
    USE comconst, ONLY : daysec
    USE comgeom, ONLY : unsaire
    USE ener, ONLY : gtot

    ! Arguments:
    REAL, intent(in):: vcov(ip1jm, llm), ucov(ip1jmp1, llm)
    real, intent(in):: teta(ip1jmp1, llm)
    REAL, intent(in):: massebx(ip1jmp1, llm), masseby(ip1jm, llm)
    real, INTENT (IN):: w(ip1jmp1, llm)
    REAL, intent(inout):: dv(ip1jm, llm), du(ip1jmp1, llm), dteta(ip1jmp1, llm)
    LOGICAL, INTENT (IN):: conser

    ! Local:
    REAL uav(ip1jmp1, llm), vav(ip1jm, llm), wsur2(ip1jmp1)
    REAL unsaire2(ip1jmp1)
    REAL deuxjour, ww, uu, vv
    INTEGER ij, l

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

    ! 2. Calculs preliminaires :

    IF (conser) THEN
       deuxjour = 2. * daysec
       unsaire2 = unsaire**2
    END IF

    ! Calcul de \bar{u}^{yy}
    DO l = 1, llm
       DO ij = iip2, ip1jmp1
          uav(ij, l) = 0.25*(ucov(ij, l)+ucov(ij-iip1, l))
       END DO
       DO ij = iip2, ip1jm
          uav(ij, l) = uav(ij, l) + uav(ij+iip1, l)
       END DO
       DO ij = 1, iip1
          uav(ij, l) = 0.
          uav(ip1jm+ij, l) = 0.
       END DO
    END DO

    ! Calcul de \bar{v}^{xx}
    DO l = 1, llm
       DO ij = 2, ip1jm
          vav(ij, l) = 0.25*(vcov(ij, l)+vcov(ij-1, l))
       END DO
       DO ij = 1, ip1jm, iip1
          vav(ij, l) = vav(ij+iim, l)
       END DO
       DO ij = 1, ip1jm - 1
          vav(ij, l) = vav(ij, l) + vav(ij+1, l)
       END DO
       DO ij = 1, ip1jm, iip1
          vav(ij+iim, l) = vav(ij, l)
       END DO
    END DO

    DO l = 1, llm - 1
       ! calcul de - w/2 au niveau l+1
       DO ij = 1, ip1jmp1
          wsur2(ij) = -0.5 * w(ij, l+1)
       END DO

       ! calcul pour "du"
       DO ij = iip2, ip1jm - 1
          ww = wsur2(ij) + wsur2(ij+1)
          uu = 0.5*(ucov(ij, l)+ucov(ij, l+1))
          du(ij, l) = du(ij, l) - ww*(uu-uav(ij, l))/massebx(ij, l)
          du(ij, l+1) = du(ij, l+1) + ww*(uu-uav(ij, l+1))/massebx(ij, l+1)
       END DO

       ! correction pour du(iip1, j, l)
       ! du(iip1, j, l)= du(1, j, l)

       DO ij = iip1 + iip1, ip1jm, iip1
          du(ij, l) = du(ij-iim, l)
          du(ij, l+1) = du(ij-iim, l+1)
       END DO

       ! calcul pour dv
       DO ij = 1, ip1jm
          ww = wsur2(ij+iip1) + wsur2(ij)
          vv = 0.5*(vcov(ij, l)+vcov(ij, l+1))
          dv(ij, l) = dv(ij, l) - ww*(vv-vav(ij, l))/masseby(ij, l)
          dv(ij, l+1) = dv(ij, l+1) + ww*(vv-vav(ij, l+1))/masseby(ij, l+1)
       END DO

       ! calcul pour dh
       ! calcul de - d(\bar{teta}^z * w) qu'on ajoute à dh
       DO ij = 1, ip1jmp1
          ww = wsur2(ij) * (teta(ij, l) + teta(ij, l + 1))
          dteta(ij, l) = dteta(ij, l) - ww
          dteta(ij, l + 1) = dteta(ij, l + 1) + ww
       end DO

       IF (conser) THEN
          gtot(l) = deuxjour * sqrt(sum(wsur2**2 * unsaire2) / ip1jmp1)
       END IF
    END DO

  END SUBROUTINE advect

end module advect_m
1	guez	46	module advect_m
2	guez	3
3	guez	25	IMPLICIT NONE
4	guez	3
5	guez	46	contains
6	guez	3
7	guez	46	SUBROUTINE advect(ucov, vcov, teta, w, massebx, masseby, du, dv, dteta, &
8			conser)
9	guez	3
10	guez	46	! From dyn3d/advect.F, version 1.1.1.1 2004/05/19 12:53:06
11			! Authors: P. Le Van , F. Hourdin
12			! Objet : calcul des termes d'advection verticale pour u, v, teta.
13			! Ces termes sont ajoutés à du, dv, dteta.
14	guez	3
15	guez	46	USE dimens_m, ONLY : iim, llm
16			USE paramet_m, ONLY : iip1, iip2, ip1jm, ip1jmp1
17			USE comconst, ONLY : daysec
18			USE comgeom, ONLY : unsaire
19			USE ener, ONLY : gtot
20	guez	3
21	guez	46	! Arguments:
22			REAL, intent(in):: vcov(ip1jm, llm), ucov(ip1jmp1, llm)
23			real, intent(in):: teta(ip1jmp1, llm)
24			REAL, intent(in):: massebx(ip1jmp1, llm), masseby(ip1jm, llm)
25			real, INTENT (IN):: w(ip1jmp1, llm)
26			REAL, intent(inout):: dv(ip1jm, llm), du(ip1jmp1, llm), dteta(ip1jmp1, llm)
27			LOGICAL, INTENT (IN):: conser
28	guez	3
29	guez	46	! Local:
30			REAL uav(ip1jmp1, llm), vav(ip1jm, llm), wsur2(ip1jmp1)
31			REAL unsaire2(ip1jmp1)
32			REAL deuxjour, ww, uu, vv
33			INTEGER ij, l
34	guez	3
35	guez	46	!-----------------------------------------------------------------------
36	guez	3
37	guez	46	! 2. Calculs preliminaires :
38	guez	3
39	guez	46	IF (conser) THEN
40			deuxjour = 2. * daysec
41			unsaire2 = unsaire**2
42			END IF
43	guez	3
44	guez	46	! Calcul de \bar{u}^{yy}
45			DO l = 1, llm
46			DO ij = iip2, ip1jmp1
47			uav(ij, l) = 0.25*(ucov(ij, l)+ucov(ij-iip1, l))
48			END DO
49			DO ij = iip2, ip1jm
50			uav(ij, l) = uav(ij, l) + uav(ij+iip1, l)
51			END DO
52			DO ij = 1, iip1
53			uav(ij, l) = 0.
54			uav(ip1jm+ij, l) = 0.
55			END DO
56			END DO
57	guez	3
58	guez	46	! Calcul de \bar{v}^{xx}
59			DO l = 1, llm
60			DO ij = 2, ip1jm
61			vav(ij, l) = 0.25*(vcov(ij, l)+vcov(ij-1, l))
62			END DO
63			DO ij = 1, ip1jm, iip1
64			vav(ij, l) = vav(ij+iim, l)
65			END DO
66			DO ij = 1, ip1jm - 1
67			vav(ij, l) = vav(ij, l) + vav(ij+1, l)
68			END DO
69			DO ij = 1, ip1jm, iip1
70			vav(ij+iim, l) = vav(ij, l)
71			END DO
72			END DO
73	guez	3
74	guez	46	DO l = 1, llm - 1
75			! calcul de - w/2 au niveau l+1
76			DO ij = 1, ip1jmp1
77			wsur2(ij) = -0.5 * w(ij, l+1)
78			END DO
79	guez	3
80	guez	46	! calcul pour "du"
81			DO ij = iip2, ip1jm - 1
82			ww = wsur2(ij) + wsur2(ij+1)
83			uu = 0.5*(ucov(ij, l)+ucov(ij, l+1))
84			du(ij, l) = du(ij, l) - ww*(uu-uav(ij, l))/massebx(ij, l)
85			du(ij, l+1) = du(ij, l+1) + ww*(uu-uav(ij, l+1))/massebx(ij, l+1)
86			END DO
87	guez	3
88	guez	46	! correction pour du(iip1, j, l)
89			! du(iip1, j, l)= du(1, j, l)
90	guez	3
91	guez	46	DO ij = iip1 + iip1, ip1jm, iip1
92			du(ij, l) = du(ij-iim, l)
93			du(ij, l+1) = du(ij-iim, l+1)
94			END DO
95	guez	3
96	guez	46	! calcul pour dv
97			DO ij = 1, ip1jm
98			ww = wsur2(ij+iip1) + wsur2(ij)
99			vv = 0.5*(vcov(ij, l)+vcov(ij, l+1))
100			dv(ij, l) = dv(ij, l) - ww*(vv-vav(ij, l))/masseby(ij, l)
101			dv(ij, l+1) = dv(ij, l+1) + ww*(vv-vav(ij, l+1))/masseby(ij, l+1)
102			END DO
103	guez	3
104	guez	46	! calcul pour dh
105			! calcul de - d(\bar{teta}^z * w) qu'on ajoute à dh
106			DO ij = 1, ip1jmp1
107			ww = wsur2(ij) * (teta(ij, l) + teta(ij, l + 1))
108			dteta(ij, l) = dteta(ij, l) - ww
109			dteta(ij, l + 1) = dteta(ij, l + 1) + ww
110			end DO
111	guez	3
112	guez	46	IF (conser) THEN
113			gtot(l) = deuxjour * sqrt(sum(wsur2*2 unsaire2) / ip1jmp1)
114			END IF
115			END DO
116	guez	3
117	guez	46	END SUBROUTINE advect
118	guez	3
119	guez	46	end module advect_m