source: trunk/LMDZ.TITAN/patch_large_domains/dyn3d/iniinterp_h.F @ 3567

Last change on this file since 3567 was 868, checked in by aslmd, 12 years ago

LMDZ.GENERIC. Added patches for large domains. Added C2H2 and C2H6 which must be named 2H2 and 2H6 in gases.def

File size: 7.4 KB
Line 
1      subroutine iniinterp_h (imo,jmo,imn,jmn ,kmax,
2     &       rlonuo,rlatvo,rlonun,rlatvn,
3     &       ktotal,iik,jjk,jk,ik,intersec,airen)
4   
5      implicit none
6
7
8
9c ---------------------------------------------------------
10c Prepare l' interpolation des variables d'une grille LMDZ
11c  dans une autre grille LMDZ en conservant la quantite
12c  totale pour les variables intensives (/m2) : ex : Pression au sol
13c
14c   (Pour chaque case autour d'un point scalaire de la nouvelle
15c    grille, on calcule la surface (en m2)en intersection avec chaque
16c    case de l'ancienne grille , pour la future interpolation)
17c
18c on calcule aussi l' aire dans la nouvelle grille
19c
20c
21c   Auteur:  F.Forget 01/1995
22c   -------
23c
24c ---------------------------------------------------------
25c   Declarations:
26c ==============
27c
28c  ARGUMENTS
29c  """""""""
30c INPUT
31       integer imo, jmo ! dimensions ancienne grille
32       integer imn,jmn  ! dimensions nouvelle grille
33       integer kmax ! taille du tableau des intersections
34       real rlonuo(imo+1)     !  Latitude et
35       real rlatvo(jmo)       !  longitude des
36       real rlonun(imn+1)     !  bord des
37       real rlatvn(jmn)     !  cases "scalaires" (input)
38
39c OUTPUT
40       integer ktotal ! nombre totale d''intersections reperees
41       integer iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
42       real intersec(kmax)  ! surface des intersections (m2)
43       real airen (imn+1,jmn+1) ! aire dans la nouvelle grille
44
45
46       
47 
48c Autres variables
49c """"""""""""""""
50       integer i,j,ii,jj
51       integer imomx,jmomx,imnmx1,jmnmx1
52!       parameter (imomx=361,jmomx=180,imnmx1=190,jmnmx1=100)
53!       parameter (imomx=361,jmomx=180,imnmx1=360,jmnmx1=190)
54       parameter (imomx=722,jmomx=360,imnmx1=720,jmnmx1=380)
55       real a(imomx+1),b(imomx+1),c(jmomx+1),d(jmomx+1)
56       real an(imnmx1+1),bn(imnmx1+1)
57       real cn(jmnmx1+1),dn(jmnmx1+1)
58       real aa, bb,cc,dd
59       real pi
60
61       pi      = 2.*ASIN( 1. )
62
63c Test dimensions imnmx1 jmnmx1
64c""""""""""""""""""""""""""""""
65c test dimensionnement tableau airetest
66      if (imn.GT.imnmx1.OR.jmn.GT.jmnmx1) then
67        write(*,*) 'STOP pb dimensionnement'
68        write(*,*) 'il faut imn < imnmx1 et jmn < jmnmx1'
69        write(*,*) 'imn imnmx1', imn,imnmx1
70        write(*,*) 'jmn jmnmx1', jmn,jmnmx1
71        call exit(1)
72      endif
73
74      if (imo.GT.imomx.OR.jmo.GT.jmomx) then
75        write(*,*) 'STOP pb dimensionnement'
76        write(*,*) 'il faut imo < imomx et jmo < jmomx'
77        write(*,*) 'imo imomx', imo,imomx
78        write(*,*) 'jmo jmomx', jmo,jmomx
79        call exit(1)
80      endif
81
82c On repere les frontieres des cases :
83c ===================================
84c
85c Attention, on ruse avec des latitudes = 90 deg au pole.
86
87
88c  Ancienne grile
89c  """"""""""""""
90      a(1) =   - rlonuo(imo+1)
91      b(1) = rlonuo(1)
92      do i=2,imo+1
93         a(i) = rlonuo(i-1)
94         b(i) =  rlonuo(i)
95      end do
96
97      d(1) = pi/2
98      do j=1,jmo
99         c(j) = rlatvo(j)
100         d(j+1) = rlatvo(j)
101      end do
102      c(jmo+1) = -pi/2
103     
104
105c  Nouvelle grille
106c  """""""""""""""
107      an(1) =  - rlonun(imn+1)
108      bn(1) = rlonun(1)
109      do i=2,imn+1
110         an(i) = rlonun(i-1)
111         bn(i) =  rlonun(i)
112      end do
113
114      dn(1) = pi/2
115      do j=1,jmn
116         cn(j) = rlatvn(j)
117         dn(j+1) = rlatvn(j)
118      end do
119      cn(jmn+1) = -pi/2
120
121c Calcul de la surface des cases scalaires de la nouvelle grille
122c ==============================================================
123      do ii=1,imn + 1
124        do jj = 1,jmn+1
125           airen(ii,jj) = (bn(ii)-an(ii))*(sin(dn(jj))-sin(cn(jj)))
126        end do
127      end do
128
129c Calcul de la surface des intersections
130c ======================================
131
132cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
133c test dimenssion kmax (calcul de ktotal)
134c"""""""""""""""""""""""""""""""""""""""
135c Calcul de ktotal, mais ralonge beaucoup en temps => pour debug
136      write(*,*)
137      write(*,*) 'DEBUT DU TEST KMAX'
138      ktotal = 0
139      do jj = 1,jmn+1
140       do j=1, jmo+1
141          if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
142              do ii=1,imn + 1
143                do i=1, imo +1
144                    if (  ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i)))
145     &        .or. ((an(ii).lt.b(i)-2*pi).and.(bn(ii).gt.a(i)-2*pi)
146     &             .and.(b(i)-2*pi.lt.-pi) )
147     &        .or. ((an(ii).lt.b(i)+2*pi).and.(bn(ii).gt.a(i)+2*pi)
148     &             .and.(a(i)+2*pi.gt.pi) )
149     &                     )then
150                      ktotal = ktotal +1
151                     end if
152                enddo
153              enddo
154           end if
155        enddo
156      enddo
157
158      if (kmax.LT.ktotal) then
159         write(*,*)
160         write(*,*) '******** ATTENTION ********'
161         write(*,*) 'kmax =',kmax
162         write(*,*) 'ktotal =',ktotal
163         write(*,*) 'Changer la valeur de kmax dans interp_horiz.F '
164         write(*,*) 'avec kmax >= ktotal'
165         write(*,*) 'EXIT dans iniinterp_h'
166         call exit(1)
167      else
168         write(*,*) 'kmax =',kmax
169         write(*,*) 'ktotal =',ktotal
170      end if
171      write(*,*) 'FIN DU TEST KMAX'
172      write(*,*)
173cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
174
175c     boucle sur la nouvelle grille
176c     """"""""""""""""""""""""""""
177      ktotal = 0
178      do jj = 1,jmn+1
179       do j=1, jmo+1
180          if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
181              do ii=1,imn + 1
182                do i=1, imo +1
183                    if (  ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i)))
184     &        .or. ((an(ii).lt.b(i)-2*pi).and.(bn(ii).gt.a(i)-2*pi)
185     &             .and.(b(i)-2*pi.lt.-pi) )
186     &        .or. ((an(ii).lt.b(i)+2*pi).and.(bn(ii).gt.a(i)+2*pi)
187     &             .and.(a(i)+2*pi.gt.pi) )
188     &                     )then
189                      ktotal = ktotal +1
190                      iik(ktotal) =ii
191                      jjk(ktotal) =jj
192                      ik(ktotal) =i
193                      jk(ktotal) =j
194
195                      dd = min(d(j), dn(jj))
196                      cc = cn(jj)
197                      if (cc.lt. c(j))cc=c(j)
198                      if((an(ii).lt.b(i)-2*pi).and.
199     &                  (bn(ii).gt.a(i)-2*pi)) then
200                          bb = min(b(i)-2*pi,bn(ii))
201                          aa = an(ii)
202                          if (aa.lt.a(i)-2*pi) aa=a(i)-2*pi
203                      else if((an(ii).lt.b(i)+2*pi).and.
204     &                       (bn(ii).gt.a(i)+2*pi)) then
205                          bb = min(b(i)+2*pi,bn(ii))
206                          aa = an(ii)
207                          if (aa.lt.a(i)+2*pi) aa=a(i)+2*pi
208                      else
209                          bb = min(b(i),bn(ii))
210                          aa = an(ii)
211                          if (aa.lt.a(i)) aa=a(i)
212                      end if
213                      intersec(ktotal)=(bb-aa)*(sin(dd)-sin(cc))
214                     end if
215                end do
216               end do
217             end if
218         end do
219       end do       
220
221
222c     TEST  INFO
223c     DO k=1,ktotal
224c      ii = iik(k)
225c      jj = jjk(k)
226c      i = ik(k)
227c      j = jk(k)
228c      if ((ii.eq.10).and.(jj.eq.10).and.(i.eq.10).and.(j.eq.10))then
229c      if (jj.eq.2.and.(ii.eq.1))then
230c          write(*,*) '**************** jj=',jj,'ii=',ii
231c          write(*,*) 'i,j =',i,j
232c          write(*,*) 'an,bn,cn,dn', an(ii), bn(ii), cn(jj),dn(jj)
233c          write(*,*) 'a,b,c,d', a(i), b(i), c(j),d(j)
234c          write(*,*) 'intersec(k)',intersec(k)
235c          write(*,*) 'airen(ii,jj)=',airen(ii,jj)
236c      end if
237c     END DO
238
239      return
240      end
Note: See TracBrowser for help on using the repository browser.