/[MITgcm]/MITgcm/pkg/profiles/profiles_interp.F
ViewVC logotype

Annotation of /MITgcm/pkg/profiles/profiles_interp.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Fri Mar 24 22:58:25 2006 UTC (18 years, 2 months ago) by heimbach
Branch: MAIN
o package cost profiles routines to better modularize them.

1 heimbach 1.1 #include "COST_CPPOPTIONS.h"
2    
3     C o==========================================================o
4     C | subroutine profiles_interp |
5     C | o 3D interpolation of model counterparts |
6     C | for netcdf profiles data |
7     C | started: Gael Forget 15-March-2006 |
8     C o==========================================================o
9    
10     SUBROUTINE profiles_interp(
11     O traj_cur_out,
12     I lon_cur,
13     I lat_cur,
14     I type_cur,
15     I file_cur,
16     I mytime,
17     I myThid
18     & )
19    
20     implicit none
21    
22     C ==================== Global Variables ===========================
23     #include "EEPARAMS.h"
24     #include "SIZE.h"
25     #include "GRID.h"
26     #include "DYNVARS.h"
27     #include "PARAMS.h"
28     #include "cal.h"
29     #include "ecco_cost.h"
30     #include "ctrl.h"
31     #include "ctrl_dummy.h"
32     #include "optim.h"
33     #ifdef ALLOW_PROFILES_CONTRIBUTION
34     # include "profiles.h"
35     #else
36     integer NLEVELMAX
37     parameter (NLEVELMAX=1)
38     #endif
39     C ==================== Routine Variables ==========================
40     _RL mytime
41     integer mythid
42     integer type_cur,file_cur
43     _RL traj_cur_out(NLEVELMAX)
44     _RL lon_cur,lat_cur
45    
46     #ifdef ALLOW_PROFILES_CONTRIBUTION
47    
48     C ==================== Local Variables ==========================
49     _RL tab_coeffs1(2,2),tab_coeffs2(2,2),tab_coeffs3(2,2)
50     _RL ponderations(2,2),pondsSUM,distance1,distance2
51     integer i,j,k,kk,kcur,iG,jG,bi,bj
52     _RL traj_cur(nR),mask_cur(nR)
53     integer prof_i,prof_j
54     _RL lon_tmp,tmp_coeff
55     c-- == end of interface ==
56    
57     DO bi = myBxLo(myThid), myBxHi(myThid)
58     DO bj = myByLo(myThid), myByHi(myThid)
59     prof_i=-10
60     prof_j=-10
61    
62     DO j=1,sNy+1
63     DO i=1,sNx+1
64     cgf value of j, south of the data point:
65     if (type_cur.NE.4) then
66     if ((yC(I,J,bi,bj).LE.lat_cur).AND.
67     &(yC(I,J+1,bi,bj).GT.lat_cur)) then
68     prof_j=j
69     endif
70     else
71     if ((yG(I,J,bi,bj).LE.lat_cur).AND.
72     &(yG(I,J+1,bi,bj).GT.lat_cur)) then
73     prof_j=j
74     endif
75     endif
76     cgf value of i, west of the data point:
77     if (type_cur.NE.3) then
78     if (xC(i+1,j,bi,bj).LT.xC(i,j,bi,bj)) then
79     lon_tmp=2*xC(i,j,bi,bj)-xC(i-1,j,bi,bj)
80     else
81     lon_tmp=xC(i+1,j,bi,bj)
82     endif
83     if ((xC(I,J,bi,bj).LE.lon_cur).AND.
84     &(lon_tmp.GT.lon_cur)) then
85     prof_i=i
86     endif
87     else
88     if (xG(i+1,j,bi,bj).LT.xG(i,j,bi,bj)) then
89     lon_tmp=2*xG(i,j,bi,bj)-xG(i-1,j,bi,bj)
90     else
91     lon_tmp=xG(i+1,j,bi,bj)
92     endif
93     if ((xG(I,J,bi,bj).LE.lon_cur).AND.
94     &(lon_tmp.GT.lon_cur)) then
95     prof_i=i
96     endif
97     endif
98     ENDDO
99     ENDDO
100    
101    
102     if ((prof_i.NE.-10).AND.(prof_j.NE.-10)) then
103     cgf) spatial interpolation coefficients
104     c meridional direction:
105     if (type_cur.NE.4) then
106     distance1=(lat_cur-yC(prof_i,prof_j,bi,bj))
107     &/(yC(prof_i,prof_j+1,bi,bj)-yC(prof_i,prof_j,bi,bj))
108     else
109     distance1=(lat_cur-yG(prof_i,prof_j,bi,bj))
110     &/(yG(prof_i,prof_j+1,bi,bj)-yG(prof_i,prof_j,bi,bj))
111     endif
112     c zonal direction:
113     if (type_cur.NE.3) then
114     if (xC(i+1,j,bi,bj).LT.xC(i,j,bi,bj)) then
115     lon_tmp=2*xC(i,j,bi,bj)-xC(i-1,j,bi,bj)
116     else
117     lon_tmp=xC(i+1,j,bi,bj)
118     endif
119     distance2=(lon_cur-xC(prof_i,prof_j,bi,bj))
120     &/(lon_tmp-xC(prof_i,prof_j,bi,bj))
121     else
122     if (xG(i+1,j,bi,bj).LT.xG(i,j,bi,bj)) then
123     lon_tmp=2*xG(i,j,bi,bj)-xG(i-1,j,bi,bj)
124     else
125     lon_tmp=xG(i+1,j,bi,bj)
126     endif
127     distance2=(lon_cur-xG(prof_i,prof_j,bi,bj))
128     &/(lon_tmp-xG(prof_i,prof_j,bi,bj))
129     endif
130    
131     tab_coeffs2(1,1)=(1-distance2)*(1-distance1)
132     tab_coeffs2(1,2)=distance2*(1-distance1)
133     tab_coeffs2(2,1)=(1-distance2)*distance1
134     tab_coeffs2(2,2)=distance2*distance1
135    
136     cgf) mask at level k:
137     do k=1,nr
138     if (type_cur.EQ.1) then
139     tab_coeffs1(1,1)=theta(prof_i,prof_j,k,bi,bj) !SO
140     tab_coeffs1(1,2)=theta(prof_i+1,prof_j,k,bi,bj) !SE
141     tab_coeffs1(2,1)=theta(prof_i,prof_j+1,k,bi,bj) !NO
142     tab_coeffs1(2,2)=theta(prof_i+1,prof_j+1,k,bi,bj) !NZ
143     tab_coeffs3(1,1)=maskC(prof_i,prof_j,k,bi,bj) !SO
144     tab_coeffs3(1,2)=maskC(prof_i+1,prof_j,k,bi,bj) !SE
145     tab_coeffs3(2,1)=maskC(prof_i,prof_j+1,k,bi,bj) !NO
146     tab_coeffs3(2,2)=maskC(prof_i+1,prof_j+1,k,bi,bj) !NZ
147     elseif (type_cur.EQ.2) then
148     tab_coeffs1(1,1)=salt(prof_i,prof_j,k,bi,bj) !SO
149     tab_coeffs1(1,2)=salt(prof_i+1,prof_j,k,bi,bj) !SE
150     tab_coeffs1(2,1)=salt(prof_i,prof_j+1,k,bi,bj) !NO
151     tab_coeffs1(2,2)=salt(prof_i+1,prof_j+1,k,bi,bj) !NZ
152     tab_coeffs3(1,1)=maskC(prof_i,prof_j,k,bi,bj) !SO
153     tab_coeffs3(1,2)=maskC(prof_i+1,prof_j,k,bi,bj) !SE
154     tab_coeffs3(2,1)=maskC(prof_i,prof_j+1,k,bi,bj) !NO
155     tab_coeffs3(2,2)=maskC(prof_i+1,prof_j+1,k,bi,bj) !NZ
156     elseif (type_cur.EQ.3) then
157     tab_coeffs1(1,1)=uVel(prof_i,prof_j,k,bi,bj) !SO
158     tab_coeffs1(1,2)=uVel(prof_i+1,prof_j,k,bi,bj) !SE
159     tab_coeffs1(2,1)=uVel(prof_i,prof_j+1,k,bi,bj) !NO
160     tab_coeffs1(2,2)=uVel(prof_i+1,prof_j+1,k,bi,bj) !NZ
161     tab_coeffs3(1,1)=maskW(prof_i,prof_j,k,bi,bj) !SO
162     tab_coeffs3(1,2)=maskW(prof_i+1,prof_j,k,bi,bj) !SE
163     tab_coeffs3(2,1)=maskW(prof_i,prof_j+1,k,bi,bj) !NO
164     tab_coeffs3(2,2)=maskW(prof_i+1,prof_j+1,k,bi,bj) !NZ
165     elseif (type_cur.EQ.4) then
166     tab_coeffs1(1,1)=vVel(prof_i,prof_j,k,bi,bj) !SO
167     tab_coeffs1(1,2)=vVel(prof_i+1,prof_j,k,bi,bj) !SE
168     tab_coeffs1(2,1)=vVel(prof_i,prof_j+1,k,bi,bj) !NO
169     tab_coeffs1(2,2)=vVel(prof_i+1,prof_j+1,k,bi,bj) !NZ
170     tab_coeffs3(1,1)=maskS(prof_i,prof_j,k,bi,bj) !SO
171     tab_coeffs3(1,2)=maskS(prof_i+1,prof_j,k,bi,bj) !SE
172     tab_coeffs3(2,1)=maskS(prof_i,prof_j+1,k,bi,bj) !NO
173     tab_coeffs3(2,2)=maskS(prof_i+1,prof_j+1,k,bi,bj) !NZ
174     else
175     tab_coeffs1(1,1)=0.
176     tab_coeffs1(2,1)=0.
177     tab_coeffs1(1,2)=0.
178     tab_coeffs1(2,2)=0.
179     tab_coeffs3(1,1)=0.
180     tab_coeffs3(2,1)=0.
181     tab_coeffs3(1,2)=0.
182     tab_coeffs3(2,2)=0.
183     endif
184    
185     ponderations(1,1)=tab_coeffs3(1,1)*tab_coeffs2(1,1)
186     ponderations(1,2)=tab_coeffs3(1,2)*tab_coeffs2(1,2)
187     ponderations(2,1)=tab_coeffs3(2,1)*tab_coeffs2(2,1)
188     ponderations(2,2)=tab_coeffs3(2,2)*tab_coeffs2(2,2)
189     pondsSUM=ponderations(1,1)+ponderations(2,1)+ponderations(1,2)+
190     & ponderations(2,2)
191     if (pondsSUM.GT.0) then
192     tab_coeffs1(1,1)=tab_coeffs1(1,1)*ponderations(1,1)/pondsSUM
193     tab_coeffs1(1,2)=tab_coeffs1(1,2)*ponderations(1,2)/pondsSUM
194     tab_coeffs1(2,1)=tab_coeffs1(2,1)*ponderations(2,1)/pondsSUM
195     tab_coeffs1(2,2)=tab_coeffs1(2,2)*ponderations(2,2)/pondsSUM
196     traj_cur(k)=tab_coeffs1(1,1)+tab_coeffs1(2,1)+
197     & tab_coeffs1(1,2)+tab_coeffs1(2,2)
198     mask_cur(k)=1
199     else
200     traj_cur(k)=0
201     mask_cur(k)=0
202     endif
203     enddo
204    
205     else
206     do k=1,nr
207     traj_cur(k)=0
208     mask_cur(k)=0
209     enddo
210     endif
211    
212     cgf vertical interpolation:
213     do kk=1,NLEVELMAX
214     traj_cur_out(kk)=0
215     prof_mask1D_cur(kk)=0
216     enddo
217     do kk=1,profdepthno(file_cur)
218     c case 1: above first grid center=> first grid center value
219     if (prof_depth(file_cur,kk).LT.-rC(1)) then
220     traj_cur_out(kk)=traj_cur(1)
221     prof_mask1D_cur(kk)=mask_cur(1)
222     c case 2: below last grid center
223     c if in cell, last cell value
224     elseif (prof_depth(file_cur,kk).GE.-rC(nr)) then
225     if ( prof_depth(file_cur,kk) .LT.
226     & (-rC(nr)+drC(nr)/2) ) then
227     traj_cur_out(kk)=traj_cur(nr)
228     prof_mask1D_cur(kk)=mask_cur(nr)
229     endif
230     c case 3: between two grid centers
231     else
232     kcur=0
233     do k=1,nr-1
234     if ((prof_depth(file_cur,kk).GE.-rC(k)).AND.
235     & (prof_depth(file_cur,kk).LT.-rC(k+1))) then
236     kcur=k
237     endif
238     enddo
239     if (kcur.EQ.0) then
240     print*,"profiles_interp unexpected case: stop 1"
241     stop
242     endif
243     if (mask_cur(kcur+1).EQ.1.) then
244     c subcase1: 2 wet points
245     c linear interpolation
246     tmp_coeff=(prof_depth(file_cur,kk)+rC(kcur))/
247     & (-rC(kcur+1)+rC(kcur))
248     traj_cur_out(kk)=(1-tmp_coeff)*traj_cur(kcur)
249     & +tmp_coeff*traj_cur(kcur+1)
250     prof_mask1D_cur(kk)=1
251     if (mask_cur(kcur).EQ.0.) then
252     print*,"profiles_interp unexpected case: stop 2"
253     stop
254     endif
255     elseif (prof_depth(file_cur,kk).LT.-rF(kcur+1)) then
256     c subcase2: only upper is wet point
257     c if in upper cell, upper cell value
258     traj_cur_out(kk)=traj_cur(kcur)
259     prof_mask1D_cur(kk)=mask_cur(kcur)
260     endif
261     endif
262     enddo
263    
264     ENDDO
265     ENDDO
266    
267     #endif
268    
269     end
270    

  ViewVC Help
Powered by ViewVC 1.1.22