Sources/dyn3d/limy.f


! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/limy.F,v 1.1.1.1 2004/05/19
! 12:53:07 lmdzadmin Exp $

SUBROUTINE limy(s0, sy, sm, pente_max)

  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget

  ! ********************************************************************
  ! Shema  d'advection " pseudo amont " .
  ! ********************************************************************
  ! q,w sont des arguments d'entree  pour le s-pg ....
  ! dq         sont des arguments de sortie pour le s-pg ....


  ! --------------------------------------------------------------------
  USE comconst
  use comgeom, only: aire
  USE conf_gcm_m
  USE dimens_m
  USE disvert_m
  USE dynetat0_m, only: rlonv, rlonu
  USE nr_util, ONLY: pi
  USE paramet_m

  IMPLICIT NONE


  ! Arguments:
  ! ----------
  REAL pente_max
  REAL s0(ip1jmp1, llm), sy(ip1jmp1, llm), sm(ip1jmp1, llm)

  ! Local
  ! ---------

  INTEGER i, ij, l

  REAL q(ip1jmp1, llm)
  REAL airej2, airejjm, airescb(iim), airesch(iim)
  REAL dyq(ip1jmp1), dyqv(ip1jm)
  REAL adyqv(ip1jm), dyqmax(ip1jmp1)

  REAL qpns, qpsn
  LOGICAL first
  SAVE first

  REAL sinlon(iip1), sinlondlon(iip1)
  REAL coslon(iip1), coslondlon(iip1)
  SAVE sinlon, coslon, sinlondlon, coslondlon


  REAL ssum
  INTEGER ismax, ismin
  EXTERNAL ssum, convflu, ismin, ismax

  DATA first/.TRUE./

  IF (first) THEN
    PRINT *, 'SCHEMA AMONT NOUVEAU'
    first = .FALSE.
    DO i = 2, iip1
      coslon(i) = cos(rlonv(i))
      sinlon(i) = sin(rlonv(i))
      coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
      sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
    END DO
    coslon(1) = coslon(iip1)
    coslondlon(1) = coslondlon(iip1)
    sinlon(1) = sinlon(iip1)
    sinlondlon(1) = sinlondlon(iip1)
  END IF


  DO l = 1, llm

    DO ij = 1, ip1jmp1
      q(ij, l) = s0(ij, l)/sm(ij, l)
      dyq(ij) = sy(ij, l)/sm(ij, l)
    END DO

    ! --------------------------------
    ! CALCUL EN LATITUDE
    ! --------------------------------

    ! On commence par calculer la valeur du traceur moyenne sur le premier
    ! cercle
    ! de latitude autour du pole (qpns pour le pole nord et qpsn pour
    ! le pole nord) qui sera utilisee pour evaluer les pentes au pole.

    airej2 = ssum(iim, aire(iip2), 1)
    airejjm = ssum(iim, aire(ip1jm-iim), 1)
    DO i = 1, iim
      airescb(i) = aire(i+iip1)*q(i+iip1, l)
      airesch(i) = aire(i+ip1jm-iip1)*q(i+ip1jm-iip1, l)
    END DO
    qpns = ssum(iim, airescb, 1)/airej2
    qpsn = ssum(iim, airesch, 1)/airejjm

    ! calcul des pentes aux points v

    DO ij = 1, ip1jm
      dyqv(ij) = q(ij, l) - q(ij+iip1, l)
      adyqv(ij) = abs(dyqv(ij))
    END DO

    ! calcul des pentes aux points scalaires

    DO ij = iip2, ip1jm
      dyqmax(ij) = min(adyqv(ij-iip1), adyqv(ij))
      dyqmax(ij) = pente_max*dyqmax(ij)
    END DO

    ! calcul des pentes aux poles

    ! calcul des pentes limites aux poles

    ! cas ou on a un extremum au pole

    ! if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    ! &   apn=0.
    ! if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    ! &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    ! &   aps=0.

    ! limitation des pentes aux poles
    ! do ij=1,iip1
    ! dyq(ij)=apn*dyq(ij)
    ! dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
    ! enddo

    ! test
    ! do ij=1,iip1
    ! dyq(iip1+ij)=0.
    ! dyq(ip1jm+ij-iip1)=0.
    ! enddo
    ! do ij=1,ip1jmp1
    ! dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
    ! enddo

    IF (dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1))<=0.) THEN
      DO ij = 1, iip1
        dyqmax(ij) = 0.
      END DO
    ELSE
      DO ij = 1, iip1
        dyqmax(ij) = pente_max*abs(dyqv(ij))
      END DO
    END IF

    IF (dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*dyqv(ismin(iim, &
        dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)<=0.) THEN
      DO ij = ip1jm + 1, ip1jmp1
        dyqmax(ij) = 0.
      END DO
    ELSE
      DO ij = ip1jm + 1, ip1jmp1
        dyqmax(ij) = pente_max*abs(dyqv(ij-iip1))
      END DO
    END IF

    ! calcul des pentes limitees

    DO ij = 1, ip1jmp1
      IF (dyqv(ij)*dyqv(ij-iip1)>0.) THEN
        dyq(ij) = sign(min(abs(dyq(ij)),dyqmax(ij)), dyq(ij))
      ELSE
        dyq(ij) = 0.
      END IF
    END DO

    DO ij = 1, ip1jmp1
      sy(ij, l) = dyq(ij)*sm(ij, l)
    END DO

  END DO ! fin de la boucle sur les couches verticales

  RETURN
END SUBROUTINE limy
1	guez	3
2	guez	81	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/limy.F,v 1.1.1.1 2004/05/19
3			! 12:53:07 lmdzadmin Exp $
4	guez	3
5	guez	81	SUBROUTINE limy(s0, sy, sm, pente_max)
6	guez	3
7	guez	81	! Auteurs: P.Le Van, F.Hourdin, F.Forget
8	guez	3
9	guez	81	! ********************************************************************
10			! Shema d'advection " pseudo amont " .
11			! ********************************************************************
12			! q,w sont des arguments d'entree pour le s-pg ....
13			! dq sont des arguments de sortie pour le s-pg ....
14	guez	3
15
16	guez	81	! --------------------------------------------------------------------
17	guez	139	USE comconst
18			use comgeom, only: aire
19			USE conf_gcm_m
20	guez	81	USE dimens_m
21			USE disvert_m
22	guez	139	USE dynetat0_m, only: rlonv, rlonu
23	guez	81	USE nr_util, ONLY: pi
24	guez	139	USE paramet_m
25
26	guez	81	IMPLICIT NONE
27	guez	3
28
29
30	guez	81	! Arguments:
31			! ----------
32			REAL pente_max
33			REAL s0(ip1jmp1, llm), sy(ip1jmp1, llm), sm(ip1jmp1, llm)
34	guez	3
35	guez	81	! Local
36			! ---------
37	guez	3
38	guez	81	INTEGER i, ij, l
39	guez	3
40	guez	81	REAL q(ip1jmp1, llm)
41			REAL airej2, airejjm, airescb(iim), airesch(iim)
42	guez	155	REAL dyq(ip1jmp1), dyqv(ip1jm)
43	guez	81	REAL adyqv(ip1jm), dyqmax(ip1jmp1)
44	guez	3
45	guez	155	REAL qpns, qpsn
46			LOGICAL first
47	guez	81	SAVE first
48	guez	3
49	guez	81	REAL sinlon(iip1), sinlondlon(iip1)
50			REAL coslon(iip1), coslondlon(iip1)
51			SAVE sinlon, coslon, sinlondlon, coslondlon
52	guez	3
53
54	guez	81	REAL ssum
55			INTEGER ismax, ismin
56			EXTERNAL ssum, convflu, ismin, ismax
57	guez	3
58	guez	81	DATA first/.TRUE./
59	guez	3
60	guez	81	IF (first) THEN
61			PRINT *, 'SCHEMA AMONT NOUVEAU'
62			first = .FALSE.
63			DO i = 2, iip1
64			coslon(i) = cos(rlonv(i))
65			sinlon(i) = sin(rlonv(i))
66			coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
67			sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
68			END DO
69			coslon(1) = coslon(iip1)
70			coslondlon(1) = coslondlon(iip1)
71			sinlon(1) = sinlon(iip1)
72			sinlondlon(1) = sinlondlon(iip1)
73			END IF
74	guez	3
75
76
77	guez	81	DO l = 1, llm
78	guez	3
79	guez	81	DO ij = 1, ip1jmp1
80			q(ij, l) = s0(ij, l)/sm(ij, l)
81			dyq(ij) = sy(ij, l)/sm(ij, l)
82			END DO
83	guez	3
84	guez	81	! --------------------------------
85			! CALCUL EN LATITUDE
86			! --------------------------------
87	guez	3
88	guez	81	! On commence par calculer la valeur du traceur moyenne sur le premier
89			! cercle
90			! de latitude autour du pole (qpns pour le pole nord et qpsn pour
91			! le pole nord) qui sera utilisee pour evaluer les pentes au pole.
92	guez	3
93	guez	81	airej2 = ssum(iim, aire(iip2), 1)
94			airejjm = ssum(iim, aire(ip1jm-iim), 1)
95			DO i = 1, iim
96			airescb(i) = aire(i+iip1)*q(i+iip1, l)
97			airesch(i) = aire(i+ip1jm-iip1)*q(i+ip1jm-iip1, l)
98			END DO
99			qpns = ssum(iim, airescb, 1)/airej2
100			qpsn = ssum(iim, airesch, 1)/airejjm
101
102			! calcul des pentes aux points v
103
104			DO ij = 1, ip1jm
105			dyqv(ij) = q(ij, l) - q(ij+iip1, l)
106			adyqv(ij) = abs(dyqv(ij))
107			END DO
108
109			! calcul des pentes aux points scalaires
110
111			DO ij = iip2, ip1jm
112			dyqmax(ij) = min(adyqv(ij-iip1), adyqv(ij))
113			dyqmax(ij) = pente_max*dyqmax(ij)
114			END DO
115
116			! calcul des pentes aux poles
117
118			! calcul des pentes limites aux poles
119
120			! cas ou on a un extremum au pole
121
122			! if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
123			! & apn=0.
124			! if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
125			! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
126			! & aps=0.
127
128			! limitation des pentes aux poles
129			! do ij=1,iip1
130			! dyq(ij)=apn*dyq(ij)
131			! dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
132			! enddo
133
134			! test
135			! do ij=1,iip1
136			! dyq(iip1+ij)=0.
137			! dyq(ip1jm+ij-iip1)=0.
138			! enddo
139			! do ij=1,ip1jmp1
140			! dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
141			! enddo
142
143			IF (dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1))<=0.) THEN
144			DO ij = 1, iip1
145			dyqmax(ij) = 0.
146			END DO
147			ELSE
148			DO ij = 1, iip1
149			dyqmax(ij) = pente_max*abs(dyqv(ij))
150			END DO
151			END IF
152
153			IF (dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*dyqv(ismin(iim, &
154			dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)<=0.) THEN
155			DO ij = ip1jm + 1, ip1jmp1
156			dyqmax(ij) = 0.
157			END DO
158			ELSE
159			DO ij = ip1jm + 1, ip1jmp1
160			dyqmax(ij) = pente_max*abs(dyqv(ij-iip1))
161			END DO
162			END IF
163
164			! calcul des pentes limitees
165
166			DO ij = 1, ip1jmp1
167			IF (dyqv(ij)*dyqv(ij-iip1)>0.) THEN
168			dyq(ij) = sign(min(abs(dyq(ij)),dyqmax(ij)), dyq(ij))
169			ELSE
170			dyq(ij) = 0.
171			END IF
172			END DO
173
174			DO ij = 1, ip1jmp1
175			sy(ij, l) = dyq(ij)*sm(ij, l)
176			END DO
177
178			END DO ! fin de la boucle sur les couches verticales
179
180			RETURN
181			END SUBROUTINE limy