libf/dyn3d/integrd.f90

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