trunk/dyn3d/integrd.f

module integrd_m

  IMPLICIT NONE

contains

  SUBROUTINE integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
       dteta, dp, vcov, ucov, teta, q, ps, masse, finvmaold, dt, leapf)

    ! From dyn3d/integrd.F, version 1.1.1.1 2004/05/19 12:53:05
    ! Auteur: P. Le Van 
    ! Objet: incrémentation des tendances dynamiques

    USE comvert, ONLY : ap, bp
    USE comgeom, ONLY : aire, apoln, apols
    USE dimens_m, ONLY : iim, jjm, llm
    USE filtreg_m, ONLY : filtreg
    use nr_util, only: assert
    USE paramet_m, ONLY : iip1, iip2, ijp1llm, ip1jm, ip1jmp1, jjp1, llmp1

    ! Arguments: 

    REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm), teta(ip1jmp1, llm)
    REAL q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nq)
    REAL ps(ip1jmp1), masse(ip1jmp1, llm)

    REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
    REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1), massem1(ip1jmp1, llm)

    REAL dv(ip1jm, llm), du(ip1jmp1, llm)
    REAL dteta(ip1jmp1, llm), dp(ip1jmp1)
    REAL finvmaold(ip1jmp1, llm)
    LOGICAL, INTENT (IN) :: leapf
    real, intent(in):: dt

    ! Local: 

    INTEGER nq
    REAL vscr(ip1jm), uscr(ip1jmp1), hscr(ip1jmp1), pscr(ip1jmp1)
    REAL massescr(ip1jmp1, llm), finvmasse(ip1jmp1, llm)
    REAL p(ip1jmp1, llmp1)
    REAL tpn, tps, tppn(iim), tpps(iim)
    REAL qpn, qps, qppn(iim), qpps(iim)
    REAL deltap(ip1jmp1, llm)

    INTEGER l, ij, iq

    REAL ssum

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

    call assert(size(q, 1) == iim + 1, size(q, 2) == jjm + 1, &
         size(q, 3) == llm, "integrd")
    nq = size(q, 4)

    DO l = 1, llm
       DO ij = 1, iip1
          ucov(ij, l) = 0.
          ucov(ij+ip1jm, l) = 0.
          uscr(ij) = 0.
          uscr(ij+ip1jm) = 0.
       END DO
    END DO

    ! integration de ps

    CALL scopy(ip1jmp1*llm, masse, 1, massescr, 1)

    DO ij = 1, ip1jmp1
       pscr(ij) = ps(ij)
       ps(ij) = psm1(ij) + dt*dp(ij)
    END DO

    DO ij = 1, ip1jmp1
       IF (ps(ij)<0.) THEN
          PRINT *, ' Au point ij = ', ij, ' , pression sol neg. ', ps(ij)
          STOP 'integrd'
       END IF
    END DO

    DO ij = 1, iim
       tppn(ij) = aire(ij)*ps(ij)
       tpps(ij) = aire(ij+ip1jm)*ps(ij+ip1jm)
    END DO
    tpn = ssum(iim, tppn, 1)/apoln
    tps = ssum(iim, tpps, 1)/apols
    DO ij = 1, iip1
       ps(ij) = tpn
       ps(ij+ip1jm) = tps
    END DO

    ! Calcul de la nouvelle masse d'air au dernier temps integre t+1

    forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
    CALL massdair(p, masse)

    CALL scopy(ijp1llm, masse, 1, finvmasse, 1)
    CALL filtreg(finvmasse, jjp1, llm, -2, 2, .TRUE., 1)

    ! integration de ucov, vcov, h

    DO l = 1, llm
       DO ij = iip2, ip1jm
          uscr(ij) = ucov(ij, l)
          ucov(ij, l) = ucovm1(ij, l) + dt*du(ij, l)
       END DO

       DO ij = 1, ip1jm
          vscr(ij) = vcov(ij, l)
          vcov(ij, l) = vcovm1(ij, l) + dt*dv(ij, l)
       END DO

       DO ij = 1, ip1jmp1
          hscr(ij) = teta(ij, l)
          teta(ij, l) = tetam1(ij, l)*massem1(ij, l)/masse(ij, l) + &
               dt*dteta(ij, l)/masse(ij, l)
       END DO

       ! Calcul de la valeur moyenne, unique aux poles pour teta

       DO ij = 1, iim
          tppn(ij) = aire(ij)*teta(ij, l)
          tpps(ij) = aire(ij+ip1jm)*teta(ij+ip1jm, l)
       END DO
       tpn = ssum(iim, tppn, 1)/apoln
       tps = ssum(iim, tpps, 1)/apols

       DO ij = 1, iip1
          teta(ij, l) = tpn
          teta(ij+ip1jm, l) = tps
       END DO

       IF (leapf) THEN
          CALL scopy(ip1jmp1, uscr(1), 1, ucovm1(1, l), 1)
          CALL scopy(ip1jm, vscr(1), 1, vcovm1(1, l), 1)
          CALL scopy(ip1jmp1, hscr(1), 1, tetam1(1, l), 1)
       END IF
    END DO

    DO l = 1, llm
       DO ij = 1, ip1jmp1
          deltap(ij, l) = p(ij, l) - p(ij, l+1)
       END DO
    END DO

    CALL qminimum(q, nq, deltap)

    ! Calcul de la valeur moyenne, unique aux poles pour q

    DO iq = 1, nq
       DO l = 1, llm
          DO ij = 1, iim
             qppn(ij) = aire(ij)*q(ij, 1, l, iq)
             qpps(ij) = aire(ij+ip1jm)*q(ij, jjm + 1, l, iq)
          END DO
          qpn = ssum(iim, qppn, 1)/apoln
          qps = ssum(iim, qpps, 1)/apols

          DO ij = 1, iip1
             q(ij, 1, l, iq) = qpn
             q(ij, jjm + 1, l, iq) = qps
          END DO
       END DO
    END DO

    CALL scopy(ijp1llm, finvmasse, 1, finvmaold, 1)

    ! Fin de l'integration de q

    IF (leapf) THEN
       CALL scopy(ip1jmp1, pscr, 1, psm1, 1)
       CALL scopy(ip1jmp1*llm, massescr, 1, massem1, 1)
    END IF

  END SUBROUTINE integrd

end module integrd_m
1	guez	32	module integrd_m
2	guez	3
3	guez	32	IMPLICIT NONE
4	guez	3
5	guez	32	contains
6	guez	3
7	guez	39	SUBROUTINE integrd(vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, &
8			dteta, dp, vcov, ucov, teta, q, ps, masse, finvmaold, dt, leapf)
9	guez	3
10	guez	39	! From dyn3d/integrd.F, version 1.1.1.1 2004/05/19 12:53:05
11			! Auteur: P. Le Van
12			! Objet: incrémentation des tendances dynamiques
13	guez	3
14	guez	32	USE comvert, ONLY : ap, bp
15			USE comgeom, ONLY : aire, apoln, apols
16	guez	40	USE dimens_m, ONLY : iim, jjm, llm
17	guez	32	USE filtreg_m, ONLY : filtreg
18	guez	39	use nr_util, only: assert
19			USE paramet_m, ONLY : iip1, iip2, ijp1llm, ip1jm, ip1jmp1, jjp1, llmp1
20	guez	3
21	guez	39	! Arguments:
22	guez	3
23	guez	39	REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm), teta(ip1jmp1, llm)
24	guez	40	REAL q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nq)
25	guez	39	REAL ps(ip1jmp1), masse(ip1jmp1, llm)
26	guez	3
27	guez	39	REAL vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
28			REAL tetam1(ip1jmp1, llm), psm1(ip1jmp1), massem1(ip1jmp1, llm)
29	guez	3
30	guez	39	REAL dv(ip1jm, llm), du(ip1jmp1, llm)
31			REAL dteta(ip1jmp1, llm), dp(ip1jmp1)
32			REAL finvmaold(ip1jmp1, llm)
33	guez	32	LOGICAL, INTENT (IN) :: leapf
34			real, intent(in):: dt
35	guez	3
36	guez	39	! Local:
37	guez	3
38	guez	39	INTEGER nq
39	guez	32	REAL vscr(ip1jm), uscr(ip1jmp1), hscr(ip1jmp1), pscr(ip1jmp1)
40	guez	39	REAL massescr(ip1jmp1, llm), finvmasse(ip1jmp1, llm)
41			REAL p(ip1jmp1, llmp1)
42	guez	32	REAL tpn, tps, tppn(iim), tpps(iim)
43			REAL qpn, qps, qppn(iim), qpps(iim)
44	guez	39	REAL deltap(ip1jmp1, llm)
45	guez	3
46	guez	32	INTEGER l, ij, iq
47	guez	3
48	guez	32	REAL ssum
49	guez	3
50	guez	32	!-----------------------------------------------------------------------
51	guez	3
52	guez	40	call assert(size(q, 1) == iim + 1, size(q, 2) == jjm + 1, &
53			size(q, 3) == llm, "integrd")
54			nq = size(q, 4)
55	guez	39
56	guez	32	DO l = 1, llm
57			DO ij = 1, iip1
58	guez	39	ucov(ij, l) = 0.
59			ucov(ij+ip1jm, l) = 0.
60	guez	32	uscr(ij) = 0.
61			uscr(ij+ip1jm) = 0.
62			END DO
63			END DO
64	guez	3
65	guez	39	! integration de ps
66	guez	3
67	guez	39	CALL scopy(ip1jmp1*llm, masse, 1, massescr, 1)
68	guez	3
69	guez	32	DO ij = 1, ip1jmp1
70			pscr(ij) = ps(ij)
71			ps(ij) = psm1(ij) + dt*dp(ij)
72			END DO
73	guez	3
74	guez	32	DO ij = 1, ip1jmp1
75			IF (ps(ij)<0.) THEN
76			PRINT *, ' Au point ij = ', ij, ' , pression sol neg. ', ps(ij)
77			STOP 'integrd'
78			END IF
79			END DO
80	guez	3
81	guez	32	DO ij = 1, iim
82			tppn(ij) = aire(ij)*ps(ij)
83			tpps(ij) = aire(ij+ip1jm)*ps(ij+ip1jm)
84			END DO
85	guez	39	tpn = ssum(iim, tppn, 1)/apoln
86			tps = ssum(iim, tpps, 1)/apols
87	guez	32	DO ij = 1, iip1
88			ps(ij) = tpn
89			ps(ij+ip1jm) = tps
90			END DO
91	guez	3
92	guez	39	! Calcul de la nouvelle masse d'air au dernier temps integre t+1
93	guez	3
94	guez	37	forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
95	guez	39	CALL massdair(p, masse)
96	guez	3
97	guez	39	CALL scopy(ijp1llm, masse, 1, finvmasse, 1)
98			CALL filtreg(finvmasse, jjp1, llm, -2, 2, .TRUE., 1)
99	guez	3
100	guez	39	! integration de ucov, vcov, h
101	guez	3
102	guez	39	DO l = 1, llm
103	guez	32	DO ij = iip2, ip1jm
104	guez	39	uscr(ij) = ucov(ij, l)
105			ucov(ij, l) = ucovm1(ij, l) + dt*du(ij, l)
106	guez	32	END DO
107	guez	3
108	guez	32	DO ij = 1, ip1jm
109	guez	39	vscr(ij) = vcov(ij, l)
110			vcov(ij, l) = vcovm1(ij, l) + dt*dv(ij, l)
111	guez	32	END DO
112	guez	3
113	guez	32	DO ij = 1, ip1jmp1
114	guez	39	hscr(ij) = teta(ij, l)
115			teta(ij, l) = tetam1(ij, l)*massem1(ij, l)/masse(ij, l) + &
116			dt*dteta(ij, l)/masse(ij, l)
117	guez	32	END DO
118	guez	3
119	guez	39	! Calcul de la valeur moyenne, unique aux poles pour teta
120	guez	3
121	guez	32	DO ij = 1, iim
122	guez	39	tppn(ij) = aire(ij)*teta(ij, l)
123			tpps(ij) = aire(ij+ip1jm)*teta(ij+ip1jm, l)
124	guez	32	END DO
125	guez	39	tpn = ssum(iim, tppn, 1)/apoln
126			tps = ssum(iim, tpps, 1)/apols
127	guez	3
128	guez	32	DO ij = 1, iip1
129	guez	39	teta(ij, l) = tpn
130			teta(ij+ip1jm, l) = tps
131	guez	32	END DO
132	guez	3
133	guez	32	IF (leapf) THEN
134	guez	39	CALL scopy(ip1jmp1, uscr(1), 1, ucovm1(1, l), 1)
135			CALL scopy(ip1jm, vscr(1), 1, vcovm1(1, l), 1)
136			CALL scopy(ip1jmp1, hscr(1), 1, tetam1(1, l), 1)
137	guez	32	END IF
138			END DO
139	guez	3
140	guez	32	DO l = 1, llm
141			DO ij = 1, ip1jmp1
142	guez	39	deltap(ij, l) = p(ij, l) - p(ij, l+1)
143	guez	32	END DO
144			END DO
145	guez	3
146	guez	39	CALL qminimum(q, nq, deltap)
147	guez	3
148	guez	39	! Calcul de la valeur moyenne, unique aux poles pour q
149	guez	3
150	guez	32	DO iq = 1, nq
151			DO l = 1, llm
152			DO ij = 1, iim
153	guez	40	qppn(ij) = aire(ij)*q(ij, 1, l, iq)
154			qpps(ij) = aire(ij+ip1jm)*q(ij, jjm + 1, l, iq)
155	guez	32	END DO
156	guez	39	qpn = ssum(iim, qppn, 1)/apoln
157			qps = ssum(iim, qpps, 1)/apols
158	guez	28
159	guez	32	DO ij = 1, iip1
160	guez	40	q(ij, 1, l, iq) = qpn
161			q(ij, jjm + 1, l, iq) = qps
162	guez	32	END DO
163			END DO
164			END DO
165	guez	28
166	guez	39	CALL scopy(ijp1llm, finvmasse, 1, finvmaold, 1)
167	guez	28
168	guez	39	! Fin de l'integration de q
169	guez	28
170	guez	32	IF (leapf) THEN
171	guez	39	CALL scopy(ip1jmp1, pscr, 1, psm1, 1)
172			CALL scopy(ip1jmp1*llm, massescr, 1, massem1, 1)
173	guez	32	END IF
174
175			END SUBROUTINE integrd
176
177			end module integrd_m