source: trunk/LMDZ.MARS/libf/dyn3d/iniinterp_h.F @ 937

Last change on this file since 937 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM

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       real a(imomx+1),b(imomx+1),c(jmomx+1),d(jmomx+1)
55       real an(imnmx1+1),bn(imnmx1+1)
56       real cn(jmnmx1+1),dn(jmnmx1+1)
57       real aa, bb,cc,dd
58       real pi
59
60       pi      = 2.*ASIN( 1. )
61
62c Test dimensions imnmx1 jmnmx1
63c""""""""""""""""""""""""""""""
64c test dimensionnement tableau airetest
65      if (imn.GT.imnmx1.OR.jmn.GT.jmnmx1) then
66        write(*,*) 'STOP pb dimensionnement'
67        write(*,*) 'il faut imn < imnmx1 et jmn < jmnmx1'
68        write(*,*) 'imn imnmx1', imn,imnmx1
69        write(*,*) 'jmn jmnmx1', jmn,jmnmx1
70        call exit(1)
71      endif
72
73      if (imo.GT.imomx.OR.jmo.GT.jmomx) then
74        write(*,*) 'STOP pb dimensionnement'
75        write(*,*) 'il faut imo < imomx et jmo < jmomx'
76        write(*,*) 'imo imomx', imo,imomx
77        write(*,*) 'jmo jmomx', jmo,jmomx
78        call exit(1)
79      endif
80
81c On repere les frontieres des cases :
82c ===================================
83c
84c Attention, on ruse avec des latitudes = 90 deg au pole.
85
86
87c  Ancienne grile
88c  """"""""""""""
89      a(1) =   - rlonuo(imo+1)
90      b(1) = rlonuo(1)
91      do i=2,imo+1
92         a(i) = rlonuo(i-1)
93         b(i) =  rlonuo(i)
94      end do
95
96      d(1) = pi/2
97      do j=1,jmo
98         c(j) = rlatvo(j)
99         d(j+1) = rlatvo(j)
100      end do
101      c(jmo+1) = -pi/2
102     
103
104c  Nouvelle grille
105c  """""""""""""""
106      an(1) =  - rlonun(imn+1)
107      bn(1) = rlonun(1)
108      do i=2,imn+1
109         an(i) = rlonun(i-1)
110         bn(i) =  rlonun(i)
111      end do
112
113      dn(1) = pi/2
114      do j=1,jmn
115         cn(j) = rlatvn(j)
116         dn(j+1) = rlatvn(j)
117      end do
118      cn(jmn+1) = -pi/2
119
120c Calcul de la surface des cases scalaires de la nouvelle grille
121c ==============================================================
122      do ii=1,imn + 1
123        do jj = 1,jmn+1
124           airen(ii,jj) = (bn(ii)-an(ii))*(sin(dn(jj))-sin(cn(jj)))
125        end do
126      end do
127
128c Calcul de la surface des intersections
129c ======================================
130
131cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
132c test dimenssion kmax (calcul de ktotal)
133c"""""""""""""""""""""""""""""""""""""""
134c Calcul de ktotal, mais ralonge beaucoup en temps => pour debug
135      write(*,*)
136      write(*,*) 'DEBUT DU TEST KMAX'
137      ktotal = 0
138      do jj = 1,jmn+1
139       do j=1, jmo+1
140          if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
141              do ii=1,imn + 1
142                do i=1, imo +1
143                    if (  ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i)))
144     &        .or. ((an(ii).lt.b(i)-2*pi).and.(bn(ii).gt.a(i)-2*pi)
145     &             .and.(b(i)-2*pi.lt.-pi) )
146     &        .or. ((an(ii).lt.b(i)+2*pi).and.(bn(ii).gt.a(i)+2*pi)
147     &             .and.(a(i)+2*pi.gt.pi) )
148     &                     )then
149                      ktotal = ktotal +1
150                     end if
151                enddo
152              enddo
153           end if
154        enddo
155      enddo
156
157      if (kmax.LT.ktotal) then
158         write(*,*)
159         write(*,*) '******** ATTENTION ********'
160         write(*,*) 'kmax =',kmax
161         write(*,*) 'ktotal =',ktotal
162         write(*,*) 'Changer la valeur de kmax dans interp_horiz.F '
163         write(*,*) 'avec kmax >= ktotal'
164         write(*,*) 'EXIT dans iniinterp_h'
165         call exit(1)
166      else
167         write(*,*) 'kmax =',kmax
168         write(*,*) 'ktotal =',ktotal
169      end if
170      write(*,*) 'FIN DU TEST KMAX'
171      write(*,*)
172cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
173
174c     boucle sur la nouvelle grille
175c     """"""""""""""""""""""""""""
176      ktotal = 0
177      do jj = 1,jmn+1
178       do j=1, jmo+1
179          if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
180              do ii=1,imn + 1
181                do i=1, imo +1
182                    if (  ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i)))
183     &        .or. ((an(ii).lt.b(i)-2*pi).and.(bn(ii).gt.a(i)-2*pi)
184     &             .and.(b(i)-2*pi.lt.-pi) )
185     &        .or. ((an(ii).lt.b(i)+2*pi).and.(bn(ii).gt.a(i)+2*pi)
186     &             .and.(a(i)+2*pi.gt.pi) )
187     &                     )then
188                      ktotal = ktotal +1
189                      iik(ktotal) =ii
190                      jjk(ktotal) =jj
191                      ik(ktotal) =i
192                      jk(ktotal) =j
193
194                      dd = min(d(j), dn(jj))
195                      cc = cn(jj)
196                      if (cc.lt. c(j))cc=c(j)
197                      if((an(ii).lt.b(i)-2*pi).and.
198     &                  (bn(ii).gt.a(i)-2*pi)) then
199                          bb = min(b(i)-2*pi,bn(ii))
200                          aa = an(ii)
201                          if (aa.lt.a(i)-2*pi) aa=a(i)-2*pi
202                      else if((an(ii).lt.b(i)+2*pi).and.
203     &                       (bn(ii).gt.a(i)+2*pi)) then
204                          bb = min(b(i)+2*pi,bn(ii))
205                          aa = an(ii)
206                          if (aa.lt.a(i)+2*pi) aa=a(i)+2*pi
207                      else
208                          bb = min(b(i),bn(ii))
209                          aa = an(ii)
210                          if (aa.lt.a(i)) aa=a(i)
211                      end if
212                      intersec(ktotal)=(bb-aa)*(sin(dd)-sin(cc))
213                     end if
214                end do
215               end do
216             end if
217         end do
218       end do       
219
220
221c     TEST  INFO
222c     DO k=1,ktotal
223c      ii = iik(k)
224c      jj = jjk(k)
225c      i = ik(k)
226c      j = jk(k)
227c      if ((ii.eq.10).and.(jj.eq.10).and.(i.eq.10).and.(j.eq.10))then
228c      if (jj.eq.2.and.(ii.eq.1))then
229c          write(*,*) '**************** jj=',jj,'ii=',ii
230c          write(*,*) 'i,j =',i,j
231c          write(*,*) 'an,bn,cn,dn', an(ii), bn(ii), cn(jj),dn(jj)
232c          write(*,*) 'a,b,c,d', a(i), b(i), c(j),d(j)
233c          write(*,*) 'intersec(k)',intersec(k)
234c          write(*,*) 'airen(ii,jj)=',airen(ii,jj)
235c      end if
236c     END DO
237
238      return
239      end
Note: See TracBrowser for help on using the repository browser.