source: trunk/LMDZ.PLUTO.old/libf/dyn3d/startsHD/interp_horiz.F @ 3607

Last change on this file since 3607 was 3175, checked in by emillour, 13 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.5 KB
Line 
1      subroutine interp_horiz (varo,varn,imo,jmo,imn,jmn,lm,
2     &  rlonuo,rlatvo,rlonun,rlatvn) 
3
4c===========================================================
5c  Interpolation Horizontales des variables d'une grille LMDZ
6c (des points SCALAIRES au point SCALAIRES)
7c  dans une autre grille LMDZ en conservant la quantite
8c  totale pour les variables intensives (/m2) : ex : Pression au sol
9c
10c Francois Forget (01/1995)
11c===========================================================
12
13      IMPLICIT NONE
14
15c   Declarations:
16c ==============
17c
18c  ARGUMENTS
19c  """""""""
20       
21       INTEGER imo, jmo ! dimensions ancienne grille (input)
22       INTEGER imn,jmn  ! dimensions nouvelle grille (input)
23
24       REAL rlonuo(imo+1)     !  Latitude et
25       REAL rlatvo(jmo)       !  longitude des
26       REAL rlonun(imn+1)     !  bord des
27       REAL rlatvn(jmn)     !  cases "scalaires" (input)
28
29       INTEGER lm ! dimension verticale (input)
30       REAL varo (imo+1, jmo+1,lm) ! var dans l'ancienne grille (input)
31       REAL varn (imn+1,jmn+1,lm) ! var dans la nouvelle grille (output)
32
33c Autres variables
34c """"""""""""""""
35       INTEGER imnmx2,jmnmx2
36c       parameter (imnmx2=190,jmnmx2=100)
37       parameter (imnmx2=360,jmnmx2=190)
38       REAL airetest(imnmx2+1,jmnmx2+1)
39       INTEGER ii,jj,l
40
41       REAL airen ((imnmx2+1)*(jmnmx2+1)) ! aire dans la nouvelle grille
42       REAL airentotn   ! aire totale pole nord dans la nouvelle grille
43       REAL airentots   ! aire totale pole sud dans la nouvelle grille
44c    Info sur les ktotal intersection entre les cases new/old grille
45
46c kmax: le nombre  max  d'intersections entre les 2 grilles horizontales
47c On fixe kmax a la taille de la grille des donnees martiennes (360x179)
48c + des pouiemes (cas ou une maille est a cheval sur 2 ou 4 mailles)
49c  Il y a un test dans iniinterp_h pour s'assurer que ktotal < kmax
50       INTEGER kmax, k, ktotal
51       !parameter (kmax = 360*179 + 200000)
52       parameter (kmax = 1440*800)
53c      parameter (kmax = 360*179 + 40000)
54
55       INTEGER iik(kmax), jjk(kmax),jk(kmax),ik(kmax)
56       REAL intersec(kmax)
57       REAL rrr
58       REAL totn, tots
59       integer prev_sumdim
60       save prev_sumdim
61       data prev_sumdim /0/
62       
63
64       logical firsttest, aire_ok
65       save firsttest
66       data firsttest /.true./
67       data aire_ok /.true./
68
69       integer imoS,jmoS,imnS,jmnS
70       save imoS,jmoS,imnS,jmnS
71       save ktotal,iik,jjk,jk,ik,intersec,airen
72       REAL pi
73
74       write(*,*) 'Start interpolation horiz'
75c Test dimensions imnmx2 jmnmx2
76c""""""""""""""""""""""""""""""
77c test dimensionnement tableau airetest
78      if (imn.GT.imnmx2.OR.jmn.GT.jmnmx2) then
79         write(*,*) 'STOP pb dimensionnement tableau airetest'
80         write(*,*) 'il faut imn < imnmx2 et jmn < jmnmx2'
81         write(*,*) 'imn imnmx2', imn,imnmx2
82         write(*,*) 'jmn jmnmx2', jmn,jmnmx2
83         call exit(1)
84      endif
85
86c initialisation
87c --------------
88c Si c'est le premier appel,  on prepare l'interpolation
89c en calculant pour chaque case autour d'un point scalaire de la
90c nouvelle grille, la surface  de intersection avec chaque
91c    case de l'ancienne grille.
92
93c  This must also be done if we change the dimension
94      if (imo+jmo+imn+jmn.ne.prev_sumdim) then
95          firsttest=.true.
96          prev_sumdim=imo+jmo+imn+jmn
97      end if       
98
99      if (firsttest) then
100        call iniinterp_h(imo,jmo,imn,jmn ,kmax,
101     &       rlonuo,rlatvo,rlonun,rlatvn,
102     &          ktotal,iik,jjk,jk,ik,intersec,airen)
103       imoS=imo
104       jmoS=jmo
105       imnS=imn
106       jmnS=jmn
107      else
108       if(imo.NE.imoS.OR.jmo.NE.jmoS.OR.imn.NE.imnS.OR.jmn.NE.jmnS) then
109        call iniinterp_h(imo,jmo,imn,jmn ,kmax,
110     &       rlonuo,rlatvo,rlonun,rlatvn,
111     &          ktotal,iik,jjk,jk,ik,intersec,airen)
112       imoS=imo
113       jmoS=jmo
114       imnS=imn
115       jmnS=jmn
116       end if
117      end if
118
119      do l=1,lm
120       do jj =1 , jmn+1
121        do ii=1, imn+1
122          varn(ii,jj,l) =0.
123        end do
124       end do
125      end do
126       
127c Interpolation horizontale
128c -------------------------
129c boucle sur toute les ktotal intersections entre les cases
130c de l'ancienne et la  nouvelle grille
131c
132      do k=1,ktotal
133        do l=1,lm
134c       Temporaire TB15 : bug interpolation at boundary long -180 and 180 pour grid 32x24
135c         if (iik(k).eq.33.or.iik(k).eq.1) then
136c          varn(iik(k),jjk(k),l) = varn(iik(k),jjk(k),l)
137c     &    + varo(ik(k), jk(k),l)*intersec(k)/airen(iik(k)
138c     &    +(jjk(k)-1)*(imn+1))*2
139c         else
140          varn(iik(k),jjk(k),l) = varn(iik(k),jjk(k),l)
141     &    + varo(ik(k), jk(k),l)*intersec(k)/airen(iik(k)
142     &    +(jjk(k)-1)*(imn+1))
143c         endif
144        end do
145      end do
146     
147c Une seule valeur au pole pour les variables ! :
148c -----------------------------------------------
149      DO l=1, lm
150         totn =0.
151         tots =0.
152
153
154c moyenne du champ au poles (ponderee par les aires)
155c"""""""""""""""""""""""""""""""
156         airentotn=0.
157         airentots=0.
158
159         do ii =1, imn+1
160            totn = totn + varn(ii,1,l)*airen(ii)
161            tots = tots + varn (ii,jmn+1,l)*airen(jmn*(imn+1)+ii)
162            airentotn=airentotn + airen(ii)
163            airentots=airentots + airen(jmn*(imn+1)+ii)
164         end do
165
166         do ii =1, imn+1
167            varn(ii,1,l) = totn/airentotn
168            varn(ii,jmn+1,l) = tots/airentots
169         end do
170
171      ENDDO
172           
173
174c---------------------------------------------------------------
175c  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST
176      if (firsttest) then
177      pi=2.*asin(1.)
178      firsttest = .false.
179      write (*,*) 'INTERP. HORIZ. : TEST SUR LES AIRES:'
180
181      do jj =1 , jmn+1
182        do ii=1, imn+1
183          airetest(ii,jj) =0.
184        end do
185      end do
186      do k=1,ktotal
187         airetest(iik(k),jjk(k))= airetest(iik(k),jjk(k)) +intersec(k)
188      end do
189      do jj =1 , jmn+1
190       do ii=1, imn+1
191         rrr = airen(ii+(jj-1)*(imn+1))/airetest(ii,jj)
192         if ((rrr.gt.1.001).or.(rrr.lt.0.999)) then
193             write (*,*) '********** PROBLEME D'' AIRES !!!',
194     &                   ' DANS L''INTERPOLATION HORIZONTALE'
195             write(*,*)'ii,jj,airen,airetest',
196     &          ii,jj,airen(ii+(jj-1)*(imn+1)),airetest(ii,jj)
197             aire_ok = .false.
198         end if
199       end do
200      end do
201      if (aire_ok) write(*,*) 'INTERP. HORIZ. : AIRES OK'
202      endif
203
204c FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST
205c --------------------------------------------------------------------
206
207
208        return
209        end
Note: See TracBrowser for help on using the repository browser.