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, 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(:, :, :) ! (ip1jmp1, 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) == ip1jmp1, size(q, 2) == llm, "integrd")
    nq = size(q, 3)

    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, l, iq)
             qpps(ij) = aire(ij+ip1jm)*q(ij+ip1jm, l, iq)
          END DO
          qpn = ssum(iim, qppn, 1)/apoln
          qps = ssum(iim, qpps, 1)/apols

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