source: LMDZ5/trunk/libf/dyn3d_common/limy.F @ 5072

Last change on this file since 5072 was 2603, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn logic.h into module logic_mod.F90
EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.9 KB
Line 
1c
2c $Id: limy.F 2603 2016-07-25 09:31:56Z abarral $
3c
4      SUBROUTINE limy(s0,sy,sm,pente_max)
5c
6c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
7c
8c    ********************************************************************
9c     Shema  d'advection " pseudo amont " .
10c    ********************************************************************
11c     q,w sont des arguments d'entree  pour le s-pg ....
12c     dq               sont des arguments de sortie pour le s-pg ....
13c
14c
15c   --------------------------------------------------------------------
16      USE comconst_mod, ONLY: pi
17      IMPLICIT NONE
18c
19      include "dimensions.h"
20      include "paramet.h"
21      include "comgeom.h"
22c
23c
24c   Arguments:
25c   ----------
26      real pente_max
27      real s0(ip1jmp1,llm),sy(ip1jmp1,llm),sm(ip1jmp1,llm)
28c
29c      Local
30c   ---------
31c
32      INTEGER i,ij,l
33c
34      REAL q(ip1jmp1,llm)
35      REAL airej2,airejjm,airescb(iim),airesch(iim)
36      real sigv,dyq(ip1jmp1),dyqv(ip1jm)
37      real adyqv(ip1jm),dyqmax(ip1jmp1)
38      REAL qbyv(ip1jm,llm)
39
40      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2
41      Logical extremum,first
42      save first
43
44      real convpn,convps,convmpn,convmps
45      real sinlon(iip1),sinlondlon(iip1)
46      real coslon(iip1),coslondlon(iip1)
47      save sinlon,coslon,sinlondlon,coslondlon
48c
49c
50      REAL      SSUM
51      integer ismax,ismin
52      EXTERNAL  SSUM, convflu,ismin,ismax
53      EXTERNAL filtreg
54
55      data first/.true./
56
57      if(first) then
58         print*,'SCHEMA AMONT NOUVEAU'
59         first=.false.
60         do i=2,iip1
61            coslon(i)=cos(rlonv(i))
62            sinlon(i)=sin(rlonv(i))
63            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
64            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
65         enddo
66         coslon(1)=coslon(iip1)
67         coslondlon(1)=coslondlon(iip1)
68         sinlon(1)=sinlon(iip1)
69         sinlondlon(1)=sinlondlon(iip1)
70      endif
71
72c
73
74      do l = 1, llm
75c
76         DO ij=1,ip1jmp1
77               q(ij,l) = s0(ij,l) / sm ( ij,l )
78               dyq(ij) = sy(ij,l) / sm ( ij,l )
79         ENDDO
80c
81c   --------------------------------
82c      CALCUL EN LATITUDE
83c   --------------------------------
84
85c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
86c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
87c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
88
89      airej2 = SSUM( iim, aire(iip2), 1 )
90      airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
91      DO i = 1, iim
92      airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
93      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
94      ENDDO
95      qpns   = SSUM( iim,  airescb ,1 ) / airej2
96      qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
97
98c   calcul des pentes aux points v
99
100      do ij=1,ip1jm
101         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
102         adyqv(ij)=abs(dyqv(ij))
103      ENDDO
104
105c   calcul des pentes aux points scalaires
106
107      do ij=iip2,ip1jm
108         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
109         dyqmax(ij)=pente_max*dyqmax(ij)
110      enddo
111
112c   calcul des pentes aux poles
113
114c   calcul des pentes limites aux poles
115
116c     print*,dyqv(iip1+1)
117c     appn=abs(dyq(1)/dyqv(iip1+1))
118c     print*,dyq(ip1jm+1)
119c     print*,dyqv(ip1jm-iip1+1)
120c     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
121c     do ij=2,iim
122c        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
123c        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
124c     enddo
125c     appn=min(pente_max/appn,1.)
126c     apps=min(pente_max/apps,1.)
127
128
129c   cas ou on a un extremum au pole
130
131c     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
132c    &   appn=0.
133c     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
134c    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
135c    &   apps=0.
136
137c   limitation des pentes aux poles
138c     do ij=1,iip1
139c        dyq(ij)=appn*dyq(ij)
140c        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
141c     enddo
142
143c   test
144c      do ij=1,iip1
145c         dyq(iip1+ij)=0.
146c         dyq(ip1jm+ij-iip1)=0.
147c      enddo
148c      do ij=1,ip1jmp1
149c         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
150c      enddo
151
152      if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
153     &   then
154         do ij=1,iip1
155            dyqmax(ij)=0.
156         enddo
157      else
158         do ij=1,iip1
159            dyqmax(ij)=pente_max*abs(dyqv(ij))
160         enddo
161      endif
162
163      if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
164     & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
165     &then
166         do ij=ip1jm+1,ip1jmp1
167            dyqmax(ij)=0.
168         enddo
169      else
170         do ij=ip1jm+1,ip1jmp1
171            dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
172         enddo
173      endif
174
175c   calcul des pentes limitees
176
177      do ij=1,ip1jmp1
178         if(dyqv(ij)*dyqv(ij-iip1).gt.0.) then
179            dyq(ij)=sign(min(abs(dyq(ij)),dyqmax(ij)),dyq(ij))
180         else
181            dyq(ij)=0.
182         endif
183      enddo
184
185         DO ij=1,ip1jmp1
186               sy(ij,l) = dyq(ij) * sm ( ij,l )
187        ENDDO
188
189      enddo ! fin de la boucle sur les couches verticales
190
191      RETURN
192      END
Note: See TracBrowser for help on using the repository browser.