source: trunk/MESOSCALE/LMDZ.MARS/libf_gcm/dyn3d/iniinterp_h.F @ 832

Last change on this file since 832 was 57, checked in by aslmd, 14 years ago

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