source: trunk/LMDZ.PLUTO.old/libf/dyn3d/iniinterp_h.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

File size: 6.1 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=1441,jmomx=720,imnmx1=1440,jmnmx1=730)
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      an(1)=-pi        !TB17
113      bn(imn+1)=pi     ! TB17
114
115      dn(1) = pi/2
116      do j=1,jmn
117         cn(j) = rlatvn(j)
118         dn(j+1) = rlatvn(j)
119      end do
120      cn(jmn+1) = -pi/2
121
122c Calcul de la surface des cases scalaires de la nouvelle grille
123c ==============================================================
124      do ii=1,imn + 1
125        do jj = 1,jmn+1
126           airen(ii,jj) = (bn(ii)-an(ii))*(sin(dn(jj))-sin(cn(jj)))
127        end do
128      end do
129
130c Calcul de la surface des intersections
131c ======================================
132
133cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
134
135c     boucle sur la nouvelle grille
136c     """"""""""""""""""""""""""""
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                      iik(ktotal) =ii
151                      jjk(ktotal) =jj
152                      ik(ktotal) =i
153                      jk(ktotal) =j
154
155                      dd = min(d(j), dn(jj))
156                      cc = cn(jj)
157                     
158                      if (cc.lt. c(j))cc=c(j)
159                      if((an(ii).lt.b(i)-2*pi).and.
160     &                  (bn(ii).gt.a(i)-2*pi)) then
161                          bb = min(b(i)-2*pi,bn(ii))
162                          aa = an(ii)
163                          if (aa.lt.a(i)-2*pi) aa=a(i)-2*pi
164                      else if((an(ii).lt.b(i)+2*pi).and.
165     &                       (bn(ii).gt.a(i)+2*pi)) then
166                          bb = min(b(i)+2*pi,bn(ii))
167                          aa = an(ii)
168                          if (aa.lt.a(i)+2*pi) aa=a(i)+2*pi
169                      else
170                          bb = min(b(i),bn(ii))
171                          aa = an(ii)
172                          if (aa.lt.a(i)) aa=a(i)
173                      end if
174                      intersec(ktotal)=(bb-aa)*(sin(dd)-sin(cc))
175c                     write(12,*) 'intersec(ktotal)',intersec(ktotal)
176                     end if
177                end do
178               end do
179             end if
180         end do
181       end do       
182
183
184c     TEST  INFO
185c       DO k=1,ktotal
186c        ii = iik(k)
187c        jj = jjk(k)
188c        i = ik(k)
189c        j = jk(k)
190c      if ((ii.eq.10).and.(jj.eq.10).and.(i.eq.10).and.(j.eq.10))then
191c      if (jj.eq.2.and.(ii.eq.1))then
192c          write(*,*) '**************** jj=',jj,'ii=',ii
193c          write(*,*) 'i,j =',i,j
194c          write(*,*) 'an,bn,cn,dn', an(ii), bn(ii), cn(jj),dn(jj)
195c          write(*,*) 'a,b,c,d', a(i), b(i), c(j),d(j)
196c          write(12,*) 'intersec(k)',intersec(k)
197c          write(13,*) 'airen(ii,jj)=',airen(ii,jj)
198c      end if
199c      end if
200c       END DO
201
202      return
203      end
Note: See TracBrowser for help on using the repository browser.