/[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.10 - (hide annotations) (download)
Mon Jul 20 19:37:09 2015 UTC (8 years, 10 months ago) by gforget
Branch: MAIN
Changes since 1.9: +19 -9 lines
- cost_profiles.F: add display of cost function terms to STDOUT
- profiles.h: increase NLEVELMAX to 110, add prof_namesmod
- profiles_init_fixed.F: use prof_names rather than hard-coded prof_T etc.
- profiles_interp.F: use prof_namesmod rather than type_cur directly.
- profiles_readparms.F: define prof_namesmod values associated with prof_names.

1 gforget 1.10 C $Header: /u/gcmpack/MITgcm/pkg/profiles/profiles_interp.F,v 1.9 2012/06/22 22:07:34 gforget Exp $
2 heimbach 1.3 C $Name: $
3    
4     #include "PROFILES_OPTIONS.h"
5 heimbach 1.1
6     C o==========================================================o
7     C | subroutine profiles_interp |
8     C | o 3D interpolation of model counterparts |
9     C | for netcdf profiles data |
10     C | started: Gael Forget 15-March-2006 |
11     C o==========================================================o
12    
13     SUBROUTINE profiles_interp(
14     O traj_cur_out,
15     I lon_cur,
16     I lat_cur,
17     I type_cur,
18     I file_cur,
19     I mytime,
20 gforget 1.5 I bi,
21     I bj,
22 heimbach 1.1 I myThid
23     & )
24    
25     implicit none
26    
27     C ==================== Global Variables ===========================
28     #include "EEPARAMS.h"
29     #include "SIZE.h"
30     #include "GRID.h"
31     #include "DYNVARS.h"
32     #include "PARAMS.h"
33 heimbach 1.3 #ifdef ALLOW_CAL
34 heimbach 1.1 #include "cal.h"
35 heimbach 1.3 #endif
36 heimbach 1.4 #ifdef ALLOW_PROFILES
37 heimbach 1.1 # include "profiles.h"
38     #else
39     integer NLEVELMAX
40     parameter (NLEVELMAX=1)
41     #endif
42 gforget 1.6 #ifdef ALLOW_PTRACERS
43     #include "PTRACERS_SIZE.h"
44 jmc 1.8 #include "PTRACERS_FIELDS.h"
45 gforget 1.6 #endif
46 heimbach 1.1 C ==================== Routine Variables ==========================
47     _RL mytime
48     integer mythid
49     integer type_cur,file_cur
50     _RL traj_cur_out(NLEVELMAX)
51     _RL lon_cur,lat_cur
52    
53 heimbach 1.4 #ifdef ALLOW_PROFILES
54 heimbach 1.1
55     C ==================== Local Variables ==========================
56     _RL tab_coeffs1(2,2),tab_coeffs2(2,2),tab_coeffs3(2,2)
57     _RL ponderations(2,2),pondsSUM,distance1,distance2
58 gforget 1.9 integer i,j,k,kk,kcur,bi,bj
59 heimbach 1.1 _RL traj_cur(nR),mask_cur(nR)
60     integer prof_i,prof_j
61 gforget 1.2 _RL lon_tmp1,lon_tmp2,lon_1,lon_2,lat_1,lat_2,tmp_coeff
62 gforget 1.10 character*(80) varname
63     integer JL
64    
65     c == external functions ==
66     integer ILNBLNK
67     EXTERNAL ILNBLNK
68    
69 heimbach 1.1 c-- == end of interface ==
70    
71 gforget 1.10 JL = ILNBLNK( prof_namesmod(type_cur) )
72     write(varname(1:80),'(1a)') prof_namesmod(type_cur)(1:JL)
73    
74 gforget 1.2 prof_i=-10
75     prof_j=-10
76     lon_1=-10
77     lon_2=-10
78     lat_1=-10
79     lat_2=-10
80 heimbach 1.1
81     DO j=1,sNy+1
82     DO i=1,sNx+1
83 gforget 1.2
84 heimbach 1.1 cgf value of j, south of the data point:
85 gforget 1.10 if (varname.NE.'vVel') then
86 gforget 1.2 if ((yC(i,j,bi,bj).LE.lat_cur).AND.
87     & (yC(i,j+1,bi,bj).GT.lat_cur)) then
88     prof_j=j
89     lat_1=yC(i,j,bi,bj)
90     lat_2=yC(i,j+1,bi,bj)
91     else
92     prof_j=prof_j
93     lat_1=lat_1
94     lat_2=lat_2
95     endif
96 heimbach 1.1 else
97 gforget 1.2 if ((yG(i,j,bi,bj).LE.lat_cur).AND.
98     & (yG(i,j+1,bi,bj).GT.lat_cur)) then
99     prof_j=j
100     lat_1=yG(i,j,bi,bj)
101     lat_2=yG(i,j+1,bi,bj)
102     else
103     prof_j=prof_j
104     lat_1=lat_1
105     lat_2=lat_2
106     endif
107 heimbach 1.1 endif
108 gforget 1.2
109 heimbach 1.1 cgf value of i, west of the data point:
110 gforget 1.10 if (varname.NE.'uVel') then
111 gforget 1.2 if (xC(i+1,j,bi,bj).LT.xC(1,j,bi,bj)) then
112     lon_tmp2=xC(i+1,j,bi,bj)+360
113 heimbach 1.1 else
114 gforget 1.2 lon_tmp2=xC(i+1,j,bi,bj)
115 heimbach 1.1 endif
116 gforget 1.2 if (xC(i,j,bi,bj).LT.xC(1,j,bi,bj)) then
117     lon_tmp1=xC(i,j,bi,bj)+360
118 heimbach 1.1 else
119 gforget 1.2 lon_tmp1=xC(i,j,bi,bj)
120 heimbach 1.1 endif
121     else
122 gforget 1.2 if (xG(i+1,j,bi,bj).LT.xG(1,j,bi,bj)) then
123     lon_tmp2=xG(i+1,j,bi,bj)+360
124 heimbach 1.1 else
125 gforget 1.2 lon_tmp2=xG(i+1,j,bi,bj)
126 heimbach 1.1 endif
127 gforget 1.2 if (xG(i,j,bi,bj).LT.xG(1,j,bi,bj)) then
128     lon_tmp1=xG(i,j,bi,bj)+360
129 heimbach 1.1 else
130 gforget 1.2 lon_tmp1=xG(i,j,bi,bj)
131 heimbach 1.1 endif
132 gforget 1.2 endif
133     if ((lon_tmp1.LE.lon_cur).AND.
134     &(lon_tmp2.GT.lon_cur)) then
135     prof_i=i
136     lon_1=lon_tmp1
137     lon_2=lon_tmp2
138     else
139     prof_i=prof_i
140     lon_1=lon_1
141     lon_2=lon_2
142 heimbach 1.1 endif
143    
144 gforget 1.2 ENDDO
145     ENDDO
146    
147    
148     if ((prof_i.NE.-10).AND.(prof_j.NE.-10)) then
149     cgf) spatial interpolation
150    
151     distance1=(lat_cur-lat_1)/(lat_2-lat_1)
152     distance2=(lon_cur-lon_1)/(lon_2-lon_1)
153     tab_coeffs2(1,1)=(1-distance2)*(1-distance1)
154     tab_coeffs2(1,2)=distance2*(1-distance1)
155     tab_coeffs2(2,1)=(1-distance2)*distance1
156     tab_coeffs2(2,2)=distance2*distance1
157 heimbach 1.1
158     do k=1,nr
159 gforget 1.10 if (varname.EQ.'theta') then
160 heimbach 1.1 tab_coeffs1(1,1)=theta(prof_i,prof_j,k,bi,bj) !SO
161     tab_coeffs1(1,2)=theta(prof_i+1,prof_j,k,bi,bj) !SE
162     tab_coeffs1(2,1)=theta(prof_i,prof_j+1,k,bi,bj) !NO
163     tab_coeffs1(2,2)=theta(prof_i+1,prof_j+1,k,bi,bj) !NZ
164     tab_coeffs3(1,1)=maskC(prof_i,prof_j,k,bi,bj) !SO
165     tab_coeffs3(1,2)=maskC(prof_i+1,prof_j,k,bi,bj) !SE
166     tab_coeffs3(2,1)=maskC(prof_i,prof_j+1,k,bi,bj) !NO
167     tab_coeffs3(2,2)=maskC(prof_i+1,prof_j+1,k,bi,bj) !NZ
168 gforget 1.10 elseif (varname.EQ.'salt') then
169 heimbach 1.1 tab_coeffs1(1,1)=salt(prof_i,prof_j,k,bi,bj) !SO
170     tab_coeffs1(1,2)=salt(prof_i+1,prof_j,k,bi,bj) !SE
171     tab_coeffs1(2,1)=salt(prof_i,prof_j+1,k,bi,bj) !NO
172     tab_coeffs1(2,2)=salt(prof_i+1,prof_j+1,k,bi,bj) !NZ
173     tab_coeffs3(1,1)=maskC(prof_i,prof_j,k,bi,bj) !SO
174     tab_coeffs3(1,2)=maskC(prof_i+1,prof_j,k,bi,bj) !SE
175     tab_coeffs3(2,1)=maskC(prof_i,prof_j+1,k,bi,bj) !NO
176     tab_coeffs3(2,2)=maskC(prof_i+1,prof_j+1,k,bi,bj) !NZ
177 gforget 1.10 elseif (varname.EQ.'uVel') then
178 heimbach 1.1 tab_coeffs1(1,1)=uVel(prof_i,prof_j,k,bi,bj) !SO
179     tab_coeffs1(1,2)=uVel(prof_i+1,prof_j,k,bi,bj) !SE
180     tab_coeffs1(2,1)=uVel(prof_i,prof_j+1,k,bi,bj) !NO
181     tab_coeffs1(2,2)=uVel(prof_i+1,prof_j+1,k,bi,bj) !NZ
182     tab_coeffs3(1,1)=maskW(prof_i,prof_j,k,bi,bj) !SO
183     tab_coeffs3(1,2)=maskW(prof_i+1,prof_j,k,bi,bj) !SE
184     tab_coeffs3(2,1)=maskW(prof_i,prof_j+1,k,bi,bj) !NO
185     tab_coeffs3(2,2)=maskW(prof_i+1,prof_j+1,k,bi,bj) !NZ
186 gforget 1.10 elseif (varname.EQ.'vVel') then
187 heimbach 1.1 tab_coeffs1(1,1)=vVel(prof_i,prof_j,k,bi,bj) !SO
188     tab_coeffs1(1,2)=vVel(prof_i+1,prof_j,k,bi,bj) !SE
189     tab_coeffs1(2,1)=vVel(prof_i,prof_j+1,k,bi,bj) !NO
190     tab_coeffs1(2,2)=vVel(prof_i+1,prof_j+1,k,bi,bj) !NZ
191     tab_coeffs3(1,1)=maskS(prof_i,prof_j,k,bi,bj) !SO
192     tab_coeffs3(1,2)=maskS(prof_i+1,prof_j,k,bi,bj) !SE
193     tab_coeffs3(2,1)=maskS(prof_i,prof_j+1,k,bi,bj) !NO
194     tab_coeffs3(2,2)=maskS(prof_i+1,prof_j+1,k,bi,bj) !NZ
195 gforget 1.10 elseif (varname.EQ.'pTracer') then
196 gforget 1.6 #ifdef ALLOW_PTRACERS
197     cgf if this gets used, an additional common block could be defined, containing
198     cgf the pTracer number (now 1, hard-coded), that would be read from the .nc input file
199     tab_coeffs1(1,1)=pTracer(prof_i,prof_j,k,bi,bj,1) !SO
200     tab_coeffs1(1,2)=pTracer(prof_i+1,prof_j,k,bi,bj,1) !SE
201     tab_coeffs1(2,1)=pTracer(prof_i,prof_j+1,k,bi,bj,1) !NO
202     tab_coeffs1(2,2)=pTracer(prof_i+1,prof_j+1,k,bi,bj,1) !NZ
203     #else
204     tab_coeffs1(1,1)=0 !SO
205     tab_coeffs1(1,2)=0 !SE
206     tab_coeffs1(2,1)=0 !NO
207     tab_coeffs1(2,2)=0 !NZ
208     #endif
209     tab_coeffs3(1,1)=maskC(prof_i,prof_j,k,bi,bj) !SO
210     tab_coeffs3(1,2)=maskC(prof_i+1,prof_j,k,bi,bj) !SE
211     tab_coeffs3(2,1)=maskC(prof_i,prof_j+1,k,bi,bj) !NO
212     tab_coeffs3(2,2)=maskC(prof_i+1,prof_j+1,k,bi,bj) !NZ
213 gforget 1.10 elseif (varname.EQ.'etaN') then
214 gforget 1.6 tab_coeffs1(1,1)=etan(prof_i,prof_j,bi,bj) !SO
215     tab_coeffs1(1,2)=etan(prof_i+1,prof_j,bi,bj) !SE
216     tab_coeffs1(2,1)=etan(prof_i,prof_j+1,bi,bj) !NO
217     tab_coeffs1(2,2)=etan(prof_i+1,prof_j+1,bi,bj) !NZ
218     tab_coeffs3(1,1)=maskC(prof_i,prof_j,1,bi,bj) !SO
219     tab_coeffs3(1,2)=maskC(prof_i+1,prof_j,1,bi,bj) !SE
220     tab_coeffs3(2,1)=maskC(prof_i,prof_j+1,1,bi,bj) !NO
221     tab_coeffs3(2,2)=maskC(prof_i+1,prof_j+1,1,bi,bj) !NZ
222 heimbach 1.1 else
223 gforget 1.2 tab_coeffs1(1,1)=0.
224     tab_coeffs1(2,1)=0.
225     tab_coeffs1(1,2)=0.
226     tab_coeffs1(2,2)=0.
227     tab_coeffs3(1,1)=0.
228     tab_coeffs3(2,1)=0.
229     tab_coeffs3(1,2)=0.
230     tab_coeffs3(2,2)=0.
231 heimbach 1.1 endif
232    
233 gforget 1.2 ponderations(1,1)=tab_coeffs3(1,1)*tab_coeffs2(1,1)
234     ponderations(1,2)=tab_coeffs3(1,2)*tab_coeffs2(1,2)
235     ponderations(2,1)=tab_coeffs3(2,1)*tab_coeffs2(2,1)
236     ponderations(2,2)=tab_coeffs3(2,2)*tab_coeffs2(2,2)
237     pondsSUM=ponderations(1,1)+ponderations(2,1)+ponderations(1,2)+
238     & ponderations(2,2)
239    
240     if (pondsSUM.GT.0) then
241 heimbach 1.1 tab_coeffs1(1,1)=tab_coeffs1(1,1)*ponderations(1,1)/pondsSUM
242     tab_coeffs1(1,2)=tab_coeffs1(1,2)*ponderations(1,2)/pondsSUM
243     tab_coeffs1(2,1)=tab_coeffs1(2,1)*ponderations(2,1)/pondsSUM
244     tab_coeffs1(2,2)=tab_coeffs1(2,2)*ponderations(2,2)/pondsSUM
245 gforget 1.2 traj_cur(k)=tab_coeffs1(1,1)+tab_coeffs1(2,1)+
246     & tab_coeffs1(1,2)+tab_coeffs1(2,2)
247     mask_cur(k)=1
248     else
249 heimbach 1.1 traj_cur(k)=0
250     mask_cur(k)=0
251 gforget 1.2 endif
252     enddo
253 heimbach 1.1
254     else
255 gforget 1.2 do k=1,nr
256 heimbach 1.1 traj_cur(k)=0
257     mask_cur(k)=0
258 gforget 1.2 enddo
259 heimbach 1.1 endif
260    
261     cgf vertical interpolation:
262     do kk=1,NLEVELMAX
263     traj_cur_out(kk)=0
264 gforget 1.5 prof_mask1D_cur(kk,bi,bj)=0
265 heimbach 1.1 enddo
266 gforget 1.5 do kk=1,ProfDepthNo(file_cur,bi,bj)
267 heimbach 1.1 c case 1: above first grid center=> first grid center value
268 gforget 1.5 if (prof_depth(file_cur,kk,bi,bj).LT.-rC(1)) then
269 gforget 1.2 traj_cur_out(kk)=traj_cur(1)
270 gforget 1.5 prof_mask1D_cur(kk,bi,bj)=mask_cur(1)
271 gforget 1.2 c case 2: just below last grid center=> last cell value
272 gforget 1.5 elseif (prof_depth(file_cur,kk,bi,bj).GE.-rC(nr)) then
273     if ( prof_depth(file_cur,kk,bi,bj) .LT.
274 gforget 1.2 & (-rC(nr)+drC(nr)/2) ) then
275     traj_cur_out(kk)=traj_cur(nr)
276 gforget 1.5 prof_mask1D_cur(kk,bi,bj)=mask_cur(nr)
277 gforget 1.2 endif
278 heimbach 1.1 c case 3: between two grid centers
279 gforget 1.2 else
280     kcur=0
281     do k=1,nr-1
282 gforget 1.5 if ((prof_depth(file_cur,kk,bi,bj).GE.-rC(k)).AND.
283     & (prof_depth(file_cur,kk,bi,bj).LT.-rC(k+1))) then
284 gforget 1.2 kcur=k
285     endif
286     enddo
287     if (kcur.EQ.0) then
288 gforget 1.7 WRITE(errorMessageUnit,'(A)')
289     & 'ERROR in PROFILES_INTERP: unexpected case 1'
290     STOP 'ABNORMAL END: S/R PROFILES_INTERP'
291 gforget 1.2 endif
292     if (mask_cur(kcur+1).EQ.1.) then
293     c subcase 1: 2 wet points=>linear interpolation
294 gforget 1.5 tmp_coeff=(prof_depth(file_cur,kk,bi,bj)+rC(kcur))/
295 gforget 1.2 & (-rC(kcur+1)+rC(kcur))
296     traj_cur_out(kk)=(1-tmp_coeff)*traj_cur(kcur)
297     & +tmp_coeff*traj_cur(kcur+1)
298 gforget 1.5 prof_mask1D_cur(kk,bi,bj)=1
299 gforget 1.2 if (mask_cur(kcur).EQ.0.) then
300 gforget 1.7 WRITE(errorMessageUnit,'(A)')
301     & 'ERROR in PROFILES_INTERP: unexpected case 2'
302     STOP 'ABNORMAL END: S/R PROFILES_INTERP'
303 gforget 1.2 endif
304 gforget 1.5 elseif (prof_depth(file_cur,kk,bi,bj).LT.-rF(kcur+1)) then
305 gforget 1.2 c subcase 2: only 1 wet point just above=>upper cell value
306     traj_cur_out(kk)=traj_cur(kcur)
307 gforget 1.5 prof_mask1D_cur(kk,bi,bj)=mask_cur(kcur)
308 gforget 1.2 endif
309     endif
310 heimbach 1.1 enddo
311    
312    
313     #endif
314    
315     end
316    

  ViewVC Help
Powered by ViewVC 1.1.22