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 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 dyq(ip1jmp1), dyqv(ip1jm)
  REAL adyqv(ip1jm), dyqmax(ip1jmp1)

  LOGICAL first
  SAVE first

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


  INTEGER ismax, ismin
  EXTERNAL 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.

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