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