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

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

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


Revision 1.10 - (hide annotations) (download)
Mon Feb 19 22:44:05 2007 UTC (17 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58w_post, checkpoint58x_post, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59, checkpoint58y_post
Changes since 1.9: +3 -3 lines
change suffix of binary output file, from ".equi.bin" to ".equi.data"

1 jmc 1.10 C $Header: /u/gcmpack/MITgcm/pkg/profiles/profiles_init_fixed.F,v 1.9 2006/12/01 00:00:47 heimbach 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_init_fixed |
8     C | o initialization for netcdf profiles data |
9     C | started: Gael Forget 15-March-2006 |
10     C o==========================================================o
11    
12     SUBROUTINE profiles_init_fixed( myThid )
13    
14     implicit none
15    
16     C ==================== Global Variables ===========================
17     #include "SIZE.h"
18     #include "EEPARAMS.h"
19     #include "PARAMS.h"
20     #include "GRID.h"
21     #include "DYNVARS.h"
22 heimbach 1.4 #ifdef ALLOW_CAL
23 heimbach 1.1 #include "cal.h"
24 heimbach 1.4 #endif
25 heimbach 1.3 #ifdef ALLOW_PROFILES
26 heimbach 1.1 # include "profiles.h"
27     # include "netcdf.inc"
28     #endif
29     C ==================== Routine Variables ==========================
30    
31     integer k,l,m,bi,bj,iG,jG, myThid,num_file,length_for_tile
32     integer fid, dimid, varid1, varid1a, varid1b
33     integer varid2,varid3
34     _RL tmpyymmdd(1000),tmphhmmss(1000),diffsecs
35     integer tmpdate(4),tmpdiff(4)
36     _RL tmp_lon, tmp_lon2(1000), tmp_lat2(1000)
37     integer vec_start(2), vec_count(2), profno_div1000, kk
38     character*(80) profilesfile, fnamedatanc
39     character*(80) fnameequinc, adfnameequinc
40     integer IL, err
41     logical exst
42    
43 heimbach 1.4 #ifdef ALLOW_PROFILES
44 heimbach 1.1
45     c == external functions ==
46     integer ILNBLNK
47 heimbach 1.8 character*(max_len_mbuf) msgbuf
48 heimbach 1.1
49     c-- == end of interface ==
50    
51 gforget 1.5 DO bi = myBxLo(myThid), myBxHi(myThid)
52     DO bj = myByLo(myThid), myByHi(myThid)
53    
54     profiles_curfile_buff(bi,bj)=0
55 heimbach 1.1
56     do m=1,NLEVELMAX
57     do l=1,1000
58 gforget 1.7 do k=1,NVARMAX
59 gforget 1.5 profiles_data_buff(m,l,k,bi,bj)=0
60     profiles_weight_buff(m,l,k,bi,bj)=0
61 heimbach 1.1 enddo
62     enddo
63     enddo
64    
65     c remplacer par une boucle ensuite :
66     do num_file=1,NFILESPROFMAX
67    
68     IL = ILNBLNK( profilesfiles(num_file) )
69     if (IL.NE.0) then
70 heimbach 1.8 write(profilesfile(1:80),'(1a)')
71     & profilesfiles(num_file)(1:IL)
72     write(msgbuf,'(a,X,i3,X,a)')
73     & 'Profiles num_file is ', num_file, profilesfile(1:80)
74     call print_message(
75     & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid)
76     else
77     write(profilesfile(1:80),'(1a)') ' '
78     write(msgbuf,'(a,X,i3,X,a)')
79     & 'Profiles num_file is ', num_file, ' empty '
80     call print_message(
81     & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid)
82 heimbach 1.1 endif
83    
84     IL = ILNBLNK( profilesfile )
85     if (IL.NE.0) then
86    
87     C===========================================================
88     c open data files and read the position vectors
89     C===========================================================
90    
91     write(fnamedatanc(1:80),'(2a)') profilesfile(1:IL),'.nc'
92 heimbach 1.8 write(msgbuf,'(a,X,i3,X,a)')
93     & 'Opening num_file ', num_file, fnamedatanc(1:80)
94     call print_message(
95     & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid)
96 gforget 1.5 err = NF_OPEN(fnamedatanc, 0, fiddata(num_file,bi,bj))
97 heimbach 1.1
98     c1) read the number of profiles :
99     cgf err = NF_OPEN(filename, 0, fid)
100 gforget 1.5 fid=fiddata(num_file,bi,bj)
101 heimbach 1.1 err = NF_INQ_DIMID(fid,'iPROF', dimid )
102 gforget 1.5 err = NF_INQ_DIMLEN(fid, dimid, ProfNo(num_file,bi,bj) )
103 heimbach 1.1 err = NF_INQ_DIMID(fid,'iDEPTH', dimid )
104     if (err.NE.NF_NOERR) then
105 heimbach 1.8 err = NF_INQ_DIMID(fid,'Z', dimid )
106 heimbach 1.1 endif
107 gforget 1.5 err = NF_INQ_DIMLEN(fid, dimid, ProfDepthNo(num_file,bi,bj) )
108 heimbach 1.8 write(msgbuf,'(a,X,4i9)')
109     & ' fid, num_file, ProfNo, ProfDepthNo ',
110     & fid, num_file, ProfNo(num_file,bi,bj),
111     & ProfDepthNo(num_file,bi,bj)
112     call print_message(
113     & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid)
114 heimbach 1.1
115     c2) read the dates and positions :
116     err = NF_INQ_VARID(fid,'depth', varid1a )
117 gforget 1.5 do k=1,ProfDepthNo(num_file,bi,bj)
118 heimbach 1.1 err = NF_GET_VAR1_DOUBLE(fid,varid1a,k,
119 gforget 1.5 & prof_depth(num_file,k,bi,bj))
120 heimbach 1.1 enddo
121    
122     err = NF_INQ_VARID(fid,'prof_YYYYMMDD', varid1a )
123     err = NF_INQ_VARID(fid,'prof_HHMMSS', varid1b )
124     err = NF_INQ_VARID(fid,'prof_lon', varid2 )
125     err = NF_INQ_VARID(fid,'prof_lat', varid3 )
126    
127 gforget 1.5 c DO bi = myBxLo(myThid), myBxHi(myThid)
128     c DO bj = myByLo(myThid), myByHi(myThid)
129 heimbach 1.1
130     do k=1,NOBSGLOB
131 gforget 1.5 prof_time(num_file,k,bi,bj)=-999
132     prof_lon(num_file,k,bi,bj)=-999
133     prof_lat(num_file,k,bi,bj)=-999
134     prof_ind_glob(num_file,k,bi,bj)=-999
135 heimbach 1.1 enddo
136    
137    
138     length_for_tile=0
139 gforget 1.5 profno_div1000=max(0,int(ProfNo(num_file,bi,bj)/1000))
140 heimbach 1.1
141 gforget 1.2 do kk=1,profno_div1000+1
142 heimbach 1.1
143 gforget 1.5 if (min(ProfNo(num_file,bi,bj), 1000*kk).GE.
144 heimbach 1.1 & 1+1000*(kk-1)) then
145    
146     vec_start(1)=1
147     vec_start(2)=1+1000*(kk-1)
148     vec_count(1)=1
149 gforget 1.5 vec_count(2)=min(1000,ProfNo(num_file,bi,bj)-1000*(kk-1))
150 heimbach 1.1
151     if ( (vec_count(2).LE.0).OR.(vec_count(2).GT.1000).OR.
152     & (vec_start(2).LE.0).OR.
153 gforget 1.5 & (vec_count(2)+vec_start(2)-1.GT.ProfNo(num_file,bi,bj)) )
154 heimbach 1.1 & then
155     print*,"stop 1",vec_start, vec_count
156     stop
157     endif
158    
159     err = NF_GET_VARA_DOUBLE(fid,varid1a,vec_start(2),
160     & vec_count(2), tmpyymmdd)
161     err = NF_GET_VARA_DOUBLE(fid,varid1b,vec_start(2),
162     & vec_count(2), tmphhmmss)
163     err = NF_GET_VARA_DOUBLE(fid,varid2,vec_start(2),
164     & vec_count(2), tmp_lon2)
165     err = NF_GET_VARA_DOUBLE(fid,varid3,vec_start(2),
166     & vec_count(2), tmp_lat2)
167    
168     if (err.NE.NF_NOERR) then
169     print*,"stop 2",vec_start(2),vec_count(2),
170 gforget 1.5 & kk,min(1000,ProfNo(num_file,bi,bj)-1000*(kk-1))
171 heimbach 1.1 stop
172     endif
173    
174 gforget 1.5 do k=1,min(1000,ProfNo(num_file,bi,bj)-1000*(kk-1))
175 heimbach 1.1
176     call cal_FullDate( int(tmpyymmdd(k)),int(tmphhmmss(k)),
177 gforget 1.5 & tmpdate,bi,bj,mythid )
178 heimbach 1.1 call cal_TimePassed( modelstartdate,tmpdate,tmpdiff,mythid )
179     call cal_ToSeconds (tmpdiff,diffsecs,mythid)
180     diffsecs=diffsecs+nIter0*deltaTclock
181    
182 gforget 1.2 if (xC(sNx+1,1,bi,bj).LT.xC(1,1,bi,bj)) then
183     tmp_lon=xC(sNx+1,1,bi,bj)+360
184 heimbach 1.1 else
185 gforget 1.2 tmp_lon=xC(sNx+1,1,bi,bj)
186 heimbach 1.1 endif
187     if ((xC(1,1,bi,bj).LE.tmp_lon2(k)).AND.
188 gforget 1.2 & (tmp_lon.GT.tmp_lon2(k)).AND.
189     & (yC(1,1,bi,bj).LE.tmp_lat2(k)).AND.
190     & (yC(1,sNy+1,bi,bj).GT.tmp_lat2(k))
191     & ) then
192     length_for_tile=length_for_tile+1
193 gforget 1.5 prof_time(num_file,length_for_tile,bi,bj)=diffsecs
194     prof_lon(num_file,length_for_tile,bi,bj)=tmp_lon2(k)
195     prof_lat(num_file,length_for_tile,bi,bj)=tmp_lat2(k)
196     prof_ind_glob(num_file,length_for_tile,bi,bj)=k+1000*(kk-1)
197 gforget 1.2 if (length_for_tile.GT.NOBSGLOB) then
198     print*,"too much profiles: need to increase NOBSGLOB,"
199     print*," or split the data file (less memory cost)"
200     stop
201     endif
202     elseif (xC(sNx+1,1,bi,bj).LT.xC(1,1,bi,bj)) then
203     if ((xC(1,1,bi,bj).LE.tmp_lon2(k)+360).AND.
204     & (tmp_lon.GT.tmp_lon2(k)+360).AND.
205 heimbach 1.1 & (yC(1,1,bi,bj).LE.tmp_lat2(k)).AND.
206     & (yC(1,sNy+1,bi,bj).GT.tmp_lat2(k))
207 gforget 1.2 & ) then
208     length_for_tile=length_for_tile+1
209 gforget 1.5 prof_time(num_file,length_for_tile,bi,bj)=diffsecs
210     prof_lon(num_file,length_for_tile,bi,bj)=tmp_lon2(k)+360
211     prof_lat(num_file,length_for_tile,bi,bj)=tmp_lat2(k)
212     prof_ind_glob(num_file,length_for_tile,bi,bj)=k+1000*(kk-1)
213 gforget 1.2 if (length_for_tile.GT.NOBSGLOB) then
214 heimbach 1.1 print*,"too much profiles: need to increase NOBSGLOB,"
215     print*," or split the data file (less memory cost)"
216     stop
217 gforget 1.2 endif
218     endif
219 heimbach 1.1 endif
220     enddo
221     endif
222     enddo
223    
224 gforget 1.5 ProfNo(num_file,bi,bj)=length_for_tile
225     print*,"fid dimid ProfNo(num_file,bi,bj)",fid, dimid,
226     & num_file, ProfNo(num_file,bi,bj)
227 heimbach 1.1
228 gforget 1.7 do k=1,NVARMAX
229 gforget 1.5 prof_num_var_cur(num_file,k,bi,bj)=0
230 heimbach 1.1 enddo
231 gforget 1.5 prof_num_var_tot(num_file,bi,bj)=0
232 heimbach 1.1
233     c3) detect available data types
234     err = NF_INQ_VARID(fid,'prof_T', varid1 )
235     if (err.EQ.NF_NOERR) then
236 gforget 1.5 vec_quantities(num_file,1,bi,bj)=.TRUE.
237     prof_num_var_tot(num_file,bi,bj)=
238     & prof_num_var_tot(num_file,bi,bj)+1
239     prof_num_var_cur(num_file,1,bi,bj)=
240     & prof_num_var_tot(num_file,bi,bj)
241 heimbach 1.1 else
242 gforget 1.5 vec_quantities(num_file,1,bi,bj)=.FALSE.
243 heimbach 1.1 endif
244     err = NF_INQ_VARID(fid,'prof_S', varid1 )
245     if (err.EQ.NF_NOERR) then
246 gforget 1.5 vec_quantities(num_file,2,bi,bj)=.TRUE.
247     prof_num_var_tot(num_file,bi,bj)=
248     & prof_num_var_tot(num_file,bi,bj)+1
249     prof_num_var_cur(num_file,2,bi,bj)=
250     & prof_num_var_tot(num_file,bi,bj)
251 heimbach 1.1 else
252 gforget 1.5 vec_quantities(num_file,2,bi,bj)=.FALSE.
253 heimbach 1.1 endif
254     err = NF_INQ_VARID(fid,'prof_U', varid1 )
255     if (err.EQ.NF_NOERR) then
256 gforget 1.5 vec_quantities(num_file,3,bi,bj)=.TRUE.
257     prof_num_var_tot(num_file,bi,bj)=
258     & prof_num_var_tot(num_file,bi,bj)+1
259     prof_num_var_cur(num_file,3,bi,bj)=
260     & prof_num_var_tot(num_file,bi,bj)
261 heimbach 1.1 else
262 gforget 1.5 vec_quantities(num_file,3,bi,bj)=.FALSE.
263 heimbach 1.1 endif
264     err = NF_INQ_VARID(fid,'prof_V', varid1 )
265     if (err.EQ.NF_NOERR) then
266 gforget 1.5 vec_quantities(num_file,4,bi,bj)=.TRUE.
267     prof_num_var_tot(num_file,bi,bj)=
268     & prof_num_var_tot(num_file,bi,bj)+1
269     prof_num_var_cur(num_file,4,bi,bj)=
270     & prof_num_var_tot(num_file,bi,bj)
271 heimbach 1.1 else
272 gforget 1.5 vec_quantities(num_file,4,bi,bj)=.FALSE.
273 heimbach 1.1 endif
274 gforget 1.6 err = NF_INQ_VARID(fid,'prof_ptr', varid1 )
275     if (err.EQ.NF_NOERR) then
276     vec_quantities(num_file,5,bi,bj)=.TRUE.
277     prof_num_var_tot(num_file,bi,bj)=
278     & prof_num_var_tot(num_file,bi,bj)+1
279     prof_num_var_cur(num_file,5,bi,bj)=
280     & prof_num_var_tot(num_file,bi,bj)
281     else
282     vec_quantities(num_file,5,bi,bj)=.FALSE.
283     endif
284     err = NF_INQ_VARID(fid,'prof_ssh', varid1 )
285     if (err.EQ.NF_NOERR) then
286     vec_quantities(num_file,6,bi,bj)=.TRUE.
287     prof_num_var_tot(num_file,bi,bj)=
288     & prof_num_var_tot(num_file,bi,bj)+1
289     prof_num_var_cur(num_file,6,bi,bj)=
290     & prof_num_var_tot(num_file,bi,bj)
291     else
292     vec_quantities(num_file,6,bi,bj)=.FALSE.
293     endif
294 heimbach 1.1
295    
296     C===========================================================
297     c create files for model counterparts to observations
298     C===========================================================
299    
300 gforget 1.5 if (ProfNo(num_file,bi,bj).GT.0) then
301 heimbach 1.1 iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
302     jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
303    
304     if (profilesfile_equi_type.EQ.1) then
305    
306     write(fnameequinc(1:80),'(2a,i3.3,a,i3.3,a)')
307     & profilesfile(1:IL),'.',iG,'.',jG,'.equi.nc'
308     write(adfnameequinc(1:80),'(3a,i3.3,a,i3.3,a)') 'ad',
309     & profilesfile(1:IL),'.',iG,'.',jG,'.equi.nc'
310    
311     inquire( file=fnameequinc, exist=exst )
312     if (.NOT.exst) then
313 gforget 1.5 call profiles_init_ncfile(num_file,
314     & fiddata(num_file,bi,bj),fnameequinc,
315     & fidforward(num_file,bi,bj),ProfNo(num_file,bi,bj),
316     & ProfDepthNo(num_file,bi,bj),
317     & bi,bj,myThid)
318     call profiles_init_ncfile(num_file,fiddata(num_file,bi,bj),
319     & adfnameequinc, fidadjoint(num_file,bi,bj),ProfNo(num_file,bi,bj),
320     & ProfDepthNo(num_file,bi,bj),bi,bj, myThid)
321 heimbach 1.1 else
322 gforget 1.5 err = NF_OPEN(fnameequinc,NF_WRITE,fidforward(num_file,bi,bj))
323     err = NF_OPEN(adfnameequinc,NF_WRITE,fidadjoint(num_file,bi,bj))
324 heimbach 1.1 endif
325    
326     else
327    
328     write(fnameequinc(1:80),'(2a,i3.3,a,i3.3,a)')
329 jmc 1.10 & profilesfile(1:IL),'.',iG,'.',jG,'.equi.data'
330 heimbach 1.1 write(adfnameequinc(1:80),'(3a,i3.3,a,i3.3,a)') 'ad',
331 jmc 1.10 & profilesfile(1:IL),'.',iG,'.',jG,'.equi.data'
332 heimbach 1.1
333     inquire( file=fnameequinc, exist=exst )
334     if (.NOT.exst) then
335 gforget 1.5 call profiles_init_ncfile(num_file,fiddata(num_file,bi,bj),
336     & fnameequinc,fidforward(num_file,bi,bj),
337     & ProfNo(num_file,bi,bj),ProfDepthNo(num_file,bi,bj),
338     & bi,bj,myThid)
339     call profiles_init_ncfile(num_file,fiddata(num_file,bi,bj),
340     & adfnameequinc, fidadjoint(num_file,bi,bj),ProfNo(num_file,bi,bj),
341     & ProfDepthNo(num_file,bi,bj),bi,bj, myThid)
342 heimbach 1.1 else
343 gforget 1.5 call MDSFINDUNIT( fidforward(num_file,bi,bj) , mythid )
344     open( fidforward(num_file,bi,bj),file=fnameequinc,
345 heimbach 1.1 & form ='unformatted',status='unknown', access='direct',
346 gforget 1.5 & recl= (ProfDepthNo(num_file,bi,bj)+1)*WORDLENGTH*2 )
347     call MDSFINDUNIT( fidadjoint(num_file,bi,bj) , mythid )
348     open( fidadjoint(num_file,bi,bj),file=adfnameequinc,
349 heimbach 1.1 & form ='unformatted',status='unknown', access='direct',
350 gforget 1.5 & recl= (ProfDepthNo(num_file,bi,bj)+1)*WORDLENGTH*2 )
351 heimbach 1.1 endif
352    
353     endif
354    
355     endif
356    
357 gforget 1.5 c ENDDO
358     c ENDDO
359 heimbach 1.1
360    
361     C===========================================================
362     else
363 gforget 1.5 ProfNo(num_file,bi,bj)=0
364 gforget 1.7 do k=1,NVARMAX
365 gforget 1.5 prof_num_var_cur(num_file,k,bi,bj)=0
366     vec_quantities(num_file,k,bi,bj)=.FALSE.
367 heimbach 1.1 enddo
368 gforget 1.5 prof_num_var_tot(num_file,bi,bj)=0
369 heimbach 1.1 do k=1,NOBSGLOB
370 gforget 1.5 prof_time(num_file,k,bi,bj)=-999
371     prof_lon(num_file,k,bi,bj)=-999
372     prof_lat(num_file,k,bi,bj)=-999
373     prof_ind_glob(num_file,k,bi,bj)=-999
374 heimbach 1.1 enddo
375    
376     endif !if (IL.NE.0) then
377     enddo ! do num_file=1,NFILESPROFMAX
378     C===========================================================
379    
380 gforget 1.5 ENDDO
381     ENDDO
382    
383 heimbach 1.1 #endif
384    
385     END
386    

  ViewVC Help
Powered by ViewVC 1.1.22