trunk/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 dimens_m
  USE paramet_m
  USE comconst
  USE disvert_m
  USE conf_gcm_m
  USE comgeom
  USE nr_util, ONLY: pi
  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 sigv, dyq(ip1jmp1), dyqv(ip1jm)
  REAL adyqv(ip1jm), dyqmax(ip1jmp1)
  REAL qbyv(ip1jm, llm)

  REAL qpns, qpsn, apn, aps, dyn1, dys1, dyn2, dys2
  LOGICAL extremum, first
  SAVE first

  REAL convpn, convps, convmpn, convmps
  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			USE dimens_m
18			USE paramet_m
19			USE comconst
20			USE disvert_m
21			USE conf_gcm_m
22			USE comgeom
23			USE nr_util, ONLY: pi
24			IMPLICIT NONE
25	guez	3
26
27
28	guez	81	! Arguments:
29			! ----------
30			REAL pente_max
31			REAL s0(ip1jmp1, llm), sy(ip1jmp1, llm), sm(ip1jmp1, llm)
32	guez	3
33	guez	81	! Local
34			! ---------
35	guez	3
36	guez	81	INTEGER i, ij, l
37	guez	3
38	guez	81	REAL q(ip1jmp1, llm)
39			REAL airej2, airejjm, airescb(iim), airesch(iim)
40			REAL sigv, dyq(ip1jmp1), dyqv(ip1jm)
41			REAL adyqv(ip1jm), dyqmax(ip1jmp1)
42			REAL qbyv(ip1jm, llm)
43	guez	3
44	guez	81	REAL qpns, qpsn, apn, aps, dyn1, dys1, dyn2, dys2
45			LOGICAL extremum, first
46			SAVE first
47	guez	3
48	guez	81	REAL convpn, convps, convmpn, convmps
49			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