/[MITgcm]/MITgcm/pkg/fizhi/update_ocean_exports.F
ViewVC logotype

Annotation of /MITgcm/pkg/fizhi/update_ocean_exports.F

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


Revision 1.27 - (hide annotations) (download)
Thu Jun 21 22:11:57 2007 UTC (17 years ago) by molod
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59g, checkpoint59f, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61l, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i
Changes since 1.26: +13 -5 lines
Add some generality to be able to read cs 612 x 102 boundary condition files

1 molod 1.27 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/update_ocean_exports.F,v 1.26 2007/03/13 01:38:18 molod Exp $
2 molod 1.1 C $Name: $
3    
4 molod 1.13 #include "FIZHI_OPTIONS.h"
5 molod 1.1 subroutine update_ocean_exports (myTime, myIter, myThid)
6     c----------------------------------------------------------------------
7     c Subroutine update_ocean_exports - 'Wrapper' routine to update
8     c the fields related to the ocean's surface that are needed
9 molod 1.3 c by fizhi (sst and sea ice extent).
10 molod 1.1 c
11     c Call: getsst (Return the current sst field-read dataset if needed)
12     c getsice (Return the current sea ice field-read data if needed)
13     c-----------------------------------------------------------------------
14     implicit none
15     #include "SIZE.h"
16     #include "GRID.h"
17 molod 1.2 #include "fizhi_ocean_coms.h"
18 molod 1.1 #include "EEPARAMS.h"
19 molod 1.3 #include "chronos.h"
20 molod 1.20 #ifdef ALLOW_EXCH2
21     #include "W2_EXCH2_TOPOLOGY.h"
22     #include "W2_EXCH2_PARAMS.h"
23     #endif /* ALLOW_EXCH2 */
24 molod 1.1
25 molod 1.14 integer myIter, myThid
26     _RL myTime
27 molod 1.1
28 molod 1.20 integer i, j, bi, bj, bislot, bjslot
29 molod 1.1 integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
30 molod 1.20 integer xsize, ysize
31 molod 1.6 _RL sstmin
32 molod 1.3 parameter ( sstmin = 273.16 )
33 molod 1.1
34 molod 1.20 #if defined(ALLOW_EXCH2)
35     PARAMETER ( xsize = exch2_domain_nxt * sNx )
36     PARAMETER ( ysize = exch2_domain_nyt * sNy )
37     #else
38     PARAMETER ( xsize = Nx )
39     PARAMETER ( ysize = Ny )
40     #endif
41     _RL sst1(xsize,ysize),sst2(xsize,ysize)
42     _RL sice1(xsize,ysize),sice2(xsize,ysize)
43 molod 1.23 integer nymd1sst(nSx,nSy),nymd2sst(nSx,nSy)
44     integer nymd1sice(nSx,nSy),nymd2sice(nSx,nSy)
45     integer nhms1sst(nSx,nSy),nhms2sst(nSx,nSy)
46     integer nhms1sice(nSx,nSy),nhms2sice(nSx,nSy)
47     integer sstdates(370,nSx,nSy),sicedates(370,nSx,nSy)
48     integer ssttimes(370,nSx,nSy),sicetimes(370,nSx,nSy)
49     logical first(nSx,nSy)
50     integer nSxnSy
51     parameter(nSxnSy = nSx*nSy)
52     data first/nSxnSy*.true./
53 molod 1.20
54 molod 1.23 save nymd1sst,nymd2sst,nymd1sice,nymd2sice
55     save nhms1sst,nhms2sst,nhms1sice,nhms2sice
56 molod 1.20 save sst1, sst2, sice1, sice2
57 molod 1.23 save sstdates, sicedates
58     save ssttimes, sicetimes
59 molod 1.20
60 molod 1.4 idim1 = 1-OLx
61     idim2 = sNx+OLx
62     jdim1 = 1-OLy
63     jdim2 = sNy+OLy
64     im1 = 1
65     im2 = sNx
66     jm1 = 1
67     jm2 = sNy
68    
69     C***********************************************************************
70 molod 1.3
71 molod 1.4 DO BJ = myByLo(myThid),myByHi(myThid)
72     DO BI = myBxLo(myThid),myBxHi(myThid)
73 molod 1.20 #if defined(ALLOW_EXCH2)
74     bislot = exch2_txglobalo(W2_myTileList(bi))-1
75     bjslot = exch2_tyglobalo(W2_myTileList(bi))-1
76     #else
77     bislot = myXGlobalLo-1+(bi-1)*sNx
78     bjslot = myYGlobalLo-1+(bj-1)*sNy
79     #endif
80 molod 1.4
81 molod 1.24 call getsst(ksst,sstclim,idim1,idim2,jdim1,jdim2,im1,im2,
82     . jm1,jm2,nSx,nSy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
83     . sst1,sst2,first(bi,bj),nymd1sst(bi,bj),nymd2sst(bi,bj),
84 molod 1.23 . nhms1sst(bi,bj),nhms2sst(bi,bj),sstdates(1,bi,bj),
85     . ssttimes(1,bi,bj),sst,myThid)
86 molod 1.24 call getsice(kice,siceclim,idim1,idim2,jdim1,jdim2,im1,im2,
87     . jm1,jm2,nSx,nSy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
88     . sice1,sice2,first(bi,bj),nymd1sice(bi,bj),nymd2sice(bi,bj),
89 molod 1.23 . nhms1sice(bi,bj),nhms2sice(bi,bj),sicedates(1,bi,bj),
90     . sicetimes(1,bi,bj),sice,myThid)
91 molod 1.3
92     c Check for Minimum Open-Water SST
93     c --------------------------------
94 molod 1.4 do j=jm1,jm2
95     do i=im1,im2
96     if(sice(i,j,bi,bj).eq.0.0 .and. sst(i,j,bi,bj).lt.sstmin)
97     . sst(i,j,bi,bj) = sstmin
98 molod 1.1 enddo
99     enddo
100    
101 molod 1.4 ENDDO
102     ENDDO
103 molod 1.22 _EXCH_XY_R8(sst,myThid)
104     _EXCH_XY_R8(sice,myThid)
105 molod 1.4
106 molod 1.1 return
107     end
108 molod 1.3
109 molod 1.25 subroutine getsice(iunit,clim,idim1,idim2,jdim1,jdim2,im1,im2,
110     . jm1,jm2,nSumx,nSumy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
111     . sicebc1,sicebc2,first,nymdbc1,nymdbc2,nhmsbc1,nhmsbc2,
112     . nymdbc,nhmsbc,sice,mythid)
113 molod 1.20 C***********************************************************************
114 molod 1.3 C
115     C!ROUTINE: GETSICE
116     C!DESCRIPTION: GETSICE returns the sea ice depth.
117     C! This routine is adaptable for any frequency
118     C! data upto a daily frequency.
119     C! note: for diurnal data ndmax should be increased.
120     C
121     C!INPUT PARAMETERS:
122 molod 1.4 C! iunit Unit number assigned to the sice data file
123     C! idim1 Start dimension in x-direction
124     C! idim2 End dimension in x-direction
125     C! jdim1 Start dimension in y-direction
126     C! jdim2 End dimension in y-direction
127     C! im1 Begin of x-direction span for filling sice
128     C! im2 End of x-direction span for filling sice
129     C! jm1 Begin of y-direction span for filling sice
130     C! jm2 End of y-direction span for filling sice
131 molod 1.11 C! nSumx Number of processors in x-direction (local processor)
132     C! nSumy Number of processors in y-direction (local processor)
133 molod 1.20 C! xsize Number of processors in x-direction (global)
134     C! ysize Number of processors in y-direction (global)
135 molod 1.4 C! bi Processor number in x-direction (local to processor)
136     C! bj Processor number in y-direction (local to processor)
137 molod 1.20 C! bislot Processor number in x-direction (global)
138     C! bjslot Processor number in y-direction (global)
139 molod 1.4 C! nymd YYMMDD of the current model timestep
140     C! nhms HHMMSS of the model time
141 molod 1.3 C
142     C!OUTPUT PARAMETERS:
143 molod 1.11 C! sice(idim1:idim2,jdim1:jdim2,nSumx,nSumy) Sea ice depth in meters
144 molod 1.3 C
145     C!ROUTINES CALLED:
146     C
147 molod 1.4 C! bcdata Reads the data for a given unit number
148     C! bcheader Reads the header info for a given unit number
149     C! interp_time Returns weights for linear interpolation
150 molod 1.3 C
151     C--------------------------------------------------------------------------
152    
153     implicit none
154 molod 1.11 #include "SIZE.h"
155 molod 1.4
156 molod 1.11 integer iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSumx,nSumy
157 molod 1.23 integer xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,mythid
158 molod 1.25 logical clim
159 molod 1.4
160 molod 1.20 _RL sicebc1(xsize,ysize)
161     _RL sicebc2(xsize,ysize)
162 molod 1.11 _RL sice(idim1:idim2,jdim1:jdim2,nSumx,nSumy)
163 molod 1.23 integer nhmsbc1,nhmsbc2,nymdbc1,nymdbc2
164     logical first
165 molod 1.3
166 molod 1.4 C Maximum number of dates in one year for the data
167     integer ndmax
168 molod 1.3 parameter (ndmax = 370)
169 molod 1.23 integer nymdbc(ndmax),nhmsbc(ndmax)
170 molod 1.3
171 molod 1.4 character*8 cname
172     character*80 cdscrip
173 molod 1.26 character*22 sicedata
174 molod 1.13 _RL fac1, fac2, lat0, lon0, timebc1, timebc2, timemod, undef
175 molod 1.23 logical found, error
176 molod 1.20 integer i,j,n,nn,iyear,iyearbc,nd,ndby3,imbc,jmbc
177 molod 1.23 integer ndatebc,nrec
178     integer nymdmod
179 molod 1.3
180    
181     C--------- Variable Initialization ---------------------------------
182    
183     data error /.false./
184    
185     c save header info
186 molod 1.23 save imbc,jmbc,lat0,lon0,ndatebc,undef
187 molod 1.3
188     c this only works for between 1950-2050
189     if (nymd .lt. 500101) then
190     nymdmod = 20000000 + nymd
191     else if (nymd .le. 991231) then
192     nymdmod = 19000000 + nymd
193     else
194     nymdmod = nymd
195     endif
196    
197 molod 1.26 iyear = nymdmod/10000
198     if(clim) then
199 molod 1.27 if(xsize.eq.192)sicedata='sice19232.weekly.clim'
200     if(xsize.eq.612)sicedata='sice612102.weekly.clim'
201 molod 1.26 else
202 molod 1.27 if(xsize.eq.192)
203     . WRITE(sicedata,'(A,I4)')'sice19232.weekly.y',iyear
204     if(xsize.eq.612)
205     . WRITE(sicedata,'(A,I4)')'sice612102.weekly.y',iyear
206 molod 1.26 endif
207 molod 1.16
208 molod 1.3 c initialize so that first time through they have values for the check
209 molod 1.23 c these values make the iyear .ne. iyearbc true anyways for
210 molod 1.3 c for the first time so first isnt checked below.
211    
212     if (first) then
213     nymdbc(2) = 0
214     nymdbc1 = 0
215     nymdbc2 = 0
216     nhmsbc1 = 0
217     nhmsbc2 = 0
218 molod 1.23 first = .false.
219 molod 1.3 endif
220    
221     C---------- Read in Header file ----------------------------------
222    
223     iyearbc = nymdbc(2)/10000
224    
225     if( iyear.ne.iyearbc ) then
226    
227 molod 1.23 close(iunit)
228 molod 1.15 open (iunit,file=sicedata,form='unformatted',access='direct',
229 molod 1.20 . recl=xsize*ysize*4)
230 molod 1.5 nrec = 1
231     call bcheader (iunit, ndmax, nrec,
232 molod 1.20 . cname, cdscrip, imbc, jmbc, lat0, lon0,
233 molod 1.5 . ndatebc, nymdbc, nhmsbc, undef, error)
234 molod 1.3
235     C--------- Check data for Compatibility ------------------------------
236    
237     C Check for correct data in boundary condition file
238     if (.not.error .and. cname.ne.'SICE') then
239 molod 1.4 write(6,*)'Wrong data in SICE boundary condition file => ',cname
240     error = .true.
241 molod 1.3 endif
242    
243     C Check Horizontal Resolution
244 molod 1.20 if(.not.error.and.imbc*jmbc.ne.xsize*ysize)then
245 molod 1.4 write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
246 molod 1.20 write(6,*) ' B.C. Resolution: ',imbc*jmbc
247     write(6,*) 'Model Resolution: ',xsize*ysize
248 molod 1.4 error = .true.
249 molod 1.3 endif
250    
251     C Check Year
252     iyearbc = nymdbc(2)/10000
253     if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
254 molod 1.5 write(6,*)' B.C. Year DOES NOT match REQUESTED Year!'
255     write(6,*)' B.C. Year: ', iyearbc
256     write(6,*)'Requested Year: ', iyear
257 molod 1.4 error = .true.
258 molod 1.3 endif
259    
260 molod 1.4 if (.not.error) then
261 molod 1.3 C if climatology, fill dates for data with current model year
262 molod 1.4 if (iyearbc.eq.0) then
263     write(6,*)
264     write(6,*) 'Climatological Dataset is being used.'
265     write(6,*) 'Current model year to be used to fill Header Dates'
266     do n = 2, ndatebc-1
267     nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
268     enddo
269 molod 1.3 C For the first date subtract 1 year from the current model NYMD
270 molod 1.4 n = 1
271     nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
272 molod 1.3 C For the last date add 1 year to the current model NYMD
273 molod 1.4 n = ndatebc
274     nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
275     endif
276 molod 1.3
277     C Write out header info
278 molod 1.23 _BEGIN_MASTER( myThid )
279 molod 1.4 write(6,*) ' Updated boundary condition data'
280     write(6,*) ' ---------------------------------'
281     write(6,*) ' Variable: ',cname
282     write(6,*) ' Description: ',cdscrip
283     write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
284     . ' Undefined value = ',undef
285     write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
286     write(6,*) ' Data valid at these times: '
287     ndby3 = ndatebc/3
288     do n = 1, ndby3*3,3
289     write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
290     1000 format(3(2x,i3,':',i8,2x,i8))
291     enddo
292     write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
293 molod 1.23 _END_MASTER( myThid )
294 molod 1.4 endif
295 molod 1.3
296 molod 1.4 endif
297 molod 1.3
298     C---------- Read sice data if necessary -------------------------------
299    
300 molod 1.4 found = .false.
301     nd = 2
302 molod 1.3
303     c If model time is not within the times of saved sice data
304     c from previous call to getsice then read new data
305    
306 molod 1.23 timemod = dfloat(nymdmod) + dfloat(nhms) /1000000
307     timebc1 = dfloat(nymdbc1) + dfloat(nhmsbc1)/1000000
308     timebc2 = dfloat(nymdbc2) + dfloat(nhmsbc2)/1000000
309 molod 1.4
310     if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then
311    
312     do while (.not.found .and. nd .le. ndatebc)
313 molod 1.23 timebc2 = dfloat(nymdbc(nd)) + dfloat(nhmsbc(nd))/1000000
314 molod 1.4 if (timebc2 .gt. timemod) then
315 molod 1.5 nymdbc1 = nymdbc(nd-1)
316     nymdbc2 = nymdbc(nd)
317     nhmsbc1 = nhmsbc(nd-1)
318     nhmsbc2 = nhmsbc(nd)
319 molod 1.20 call bcdata (iunit,imbc,jmbc,nd,nd+1,sicebc1,sicebc2)
320 molod 1.5 found = .true.
321 molod 1.4 else
322 molod 1.5 nd = nd + 1
323 molod 1.4 endif
324     enddo
325 molod 1.3
326     c Otherwise the data from the last time in getsice surrounds the
327     c current model time.
328    
329 molod 1.4 else
330     found = .true.
331     endif
332    
333     if (.not.found) then
334     print *, 'STOP: Could not find SICE dates for model time.'
335     call my_finalize
336     call my_exit (101)
337     endif
338 molod 1.3
339     C---------- Interpolate sice data ------------------------------------
340    
341 molod 1.4 call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
342     . fac1,fac2)
343 molod 1.3
344 molod 1.4 do j = jm1,jm2
345     do i = im1,im2
346 molod 1.20 sice(i,j,bi,bj) = sicebc1(i+bislot,j+bjslot)*fac1
347     . + sicebc2(i+bislot,j+bjslot)*fac2
348 molod 1.3 c average to 0 or 1
349     c -----------------
350 molod 1.4 if (sice(i,j,bi,bj) .ge. 0.5) then
351     sice(i,j,bi,bj) = 1.
352     else
353     sice(i,j,bi,bj) = 0.
354     endif
355     enddo
356     enddo
357 molod 1.3
358     C---------- Fill sice with depth of ice ------------------------------------
359 molod 1.4 do j = jm1,jm2
360     do i = im1,im2
361     if (sice(i,j,bi,bj) .eq. 1.) then
362     sice(i,j,bi,bj) = 3.
363     endif
364     enddo
365     enddo
366 molod 1.3 C---------------------------------------------------------------------------
367    
368     return
369     end
370 molod 1.25 subroutine getsst(iunit,clim,idim1,idim2,jdim1,jdim2,im1,im2,
371     . jm1,jm2,nSumx,nSumy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
372     . sstbc1,sstbc2,first,nymdbc1,nymdbc2,nhmsbc1,nhmsbc2,
373     . nymdbc,nhmsbc,sst,mythid)
374 molod 1.20 C***********************************************************************
375 molod 1.3 C
376     C!ROUTINE: GETSST
377     C!DESCRIPTION: GETSST gets the SST data.
378     C! This routine is adaptable for any frequency
379     C! data upto a daily frequency.
380     C! note: for diurnal data ndmax should be increased.
381     C
382     C!INPUT PARAMETERS:
383 molod 1.5 C! iunit Unit number assigned to the sice data file
384     C! idim1 Start dimension in x-direction
385     C! idim2 End dimension in x-direction
386     C! jdim1 Start dimension in y-direction
387     C! jdim2 End dimension in y-direction
388 molod 1.20 C! im1 Begin of x-direction span for filling sst
389     C! im2 End of x-direction span for filling sst
390     C! jm1 Begin of y-direction span for filling sst
391     C! jm2 End of y-direction span for filling sst
392 molod 1.11 C! nSumx Number of processors in x-direction (local processor)
393     C! nSumy Number of processors in y-direction (local processor)
394 molod 1.20 C! xsize x-dimension of global array
395     C! ysize y-dimension of global array
396 molod 1.5 C! bi Processor number in x-direction (local to processor)
397     C! bj Processor number in y-direction (local to processor)
398 molod 1.20 C! bislot Slot number into global array in x-direction (global)
399     C! bjslot Slot number into global array in y-direction (global)
400 molod 1.5 C! nymd YYMMDD of the current model timestep
401     C! nhms HHMMSS of the model time
402 molod 1.3 C
403     C!OUTPUT PARAMETERS:
404 molod 1.11 C! sst(idim1:idim2,jdim1:jdim2,nSumx,nSumy) Sea surface temperature (K)
405 molod 1.3 C
406     C!ROUTINES CALLED:
407     C
408     C! bcdata Reads the data for a given unit number
409     C! bcheader Reads the header info for a given unit number
410 molod 1.20 C! interp_time Returns weights for linear interpolation
411 molod 1.3 C
412     C--------------------------------------------------------------------------
413    
414     implicit none
415 molod 1.11 #include "SIZE.h"
416 molod 1.5
417 molod 1.11 integer iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSumx,nSumy
418 molod 1.23 integer xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,mythid
419 molod 1.25 logical clim
420 molod 1.3
421 molod 1.20 _RL sstbc1(xsize,ysize)
422     _RL sstbc2(xsize,ysize)
423 molod 1.11 _RL sst(idim1:idim2,jdim1:jdim2,nSumx,nSumy)
424 molod 1.23 integer nhmsbc1,nhmsbc2,nymdbc1,nymdbc2
425     logical first
426 molod 1.3
427 molod 1.5 C Maximum number of dates in one year for the data
428     integer ndmax
429 molod 1.3 parameter (ndmax = 370)
430 molod 1.23 integer nymdbc(ndmax),nhmsbc(ndmax)
431 molod 1.3
432 molod 1.5 character*8 cname
433     character*80 cdscrip
434 molod 1.26 character*21 sstdata
435 molod 1.13 _RL fac1, fac2, lat0, lon0, timebc1, timebc2, timemod, undef
436 molod 1.23 logical found, error
437 molod 1.20 integer i,j,n,nn,iyear,iyearbc,nd,ndby3,imbc,jmbc
438 molod 1.23 integer ndatebc,nrec
439     integer nymdmod
440 molod 1.5
441 molod 1.3
442     C--------- Variable Initialization ---------------------------------
443    
444     data error /.false./
445    
446     c save header info
447 molod 1.23 save imbc,jmbc,lat0,lon0,ndatebc,undef
448 molod 1.3
449     c this only works for between 1950-2050
450     if (nymd .lt. 500101) then
451     nymdmod = 20000000 + nymd
452     else if (nymd .le. 991231) then
453     nymdmod = 19000000 + nymd
454     else
455     nymdmod = nymd
456     endif
457    
458 molod 1.26 iyear = nymdmod/10000
459     if(clim) then
460 molod 1.27 if(xsize.eq.192)sstdata='sst19232.weekly.clim'
461     if(xsize.eq.612)sstdata='sst612102.weekly.clim'
462 molod 1.26 else
463 molod 1.27 if(xsize.eq.192)
464     . WRITE(sstdata,'(A,I4)')'sst19232.weekly.y',iyear
465     if(xsize.eq.612)
466     . WRITE(sstdata,'(A,I4)')'sst612102.weekly.y',iyear
467 molod 1.26 endif
468 molod 1.16
469 molod 1.3 c initialize so that first time through they have values for the check
470     c these vaules make the iyear .ne. iyearbc true anyways for
471     c for the first time so first isnt checked below.
472     if (first) then
473     nymdbc(2) = 0
474     nymdbc1 = 0
475     nymdbc2 = 0
476     nhmsbc1 = 0
477     nhmsbc2 = 0
478 molod 1.23 first = .false.
479 molod 1.3 endif
480    
481     C---------- Read in Header file ----------------------------------
482    
483     iyearbc = nymdbc(2)/10000
484    
485     if( iyear.ne.iyearbc ) then
486    
487 molod 1.23 close(iunit)
488 molod 1.15 open (iunit,file=sstdata,form='unformatted',access='direct',
489 molod 1.20 . recl=xsize*ysize*4)
490 molod 1.5 nrec = 1
491     call bcheader (iunit, ndmax, nrec,
492 molod 1.20 . cname, cdscrip, imbc, jmbc, lat0, lon0,
493 molod 1.5 . ndatebc, nymdbc, nhmsbc, undef, error)
494 molod 1.3
495     C--------- Check data for Compatibility
496    
497     C Check for correct data in boundary condition file
498     if (.not.error .and. cname.ne.'SST') then
499 molod 1.5 write(6,*)'Wrong data in SST boundary condition file => ',cname
500     error = .true.
501 molod 1.3 endif
502    
503     C Check Horizontal Resolution
504 molod 1.20 if(.not.error.and.imbc*jmbc.ne.xsize*ysize)then
505 molod 1.5 write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
506 molod 1.20 write(6,*) ' B.C. Resolution: ',imbc*jmbc
507     write(6,*) 'Model Resolution: ',xsize*ysize
508 molod 1.5 error = .true.
509 molod 1.3 endif
510    
511     C Check Year
512     iyearbc = nymdbc(2)/10000
513     if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
514 molod 1.5 write(6,*)' B.C. Year DOES NOT match REQUESTED Year!'
515     write(6,*)' B.C. Year: ', iyearbc
516     write(6,*)'Requested Year: ', iyear
517     error = .true.
518 molod 1.3 endif
519    
520     if (.not.error) then
521     C if climatology, fill dates for data with current model year
522 molod 1.5 if (iyearbc.eq.0) then
523     write(6,*)
524 molod 1.9 write(6,*)'Climatological Dataset is being used.'
525     write(6,*)'Current model year is used to fill Header Dates'
526 molod 1.5 do n = 2, ndatebc-1
527     nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
528     enddo
529 molod 1.3 C For the first date subtract 1 year from the current model NYMD
530 molod 1.5 n = 1
531     nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
532 molod 1.3 C For the last date add 1 year to the current model NYMD
533 molod 1.5 n = ndatebc
534     nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
535     endif
536 molod 1.3
537     C Write out header info
538 molod 1.23 _BEGIN_MASTER( myThid )
539 molod 1.5 write(6,*) ' Updated boundary condition data'
540     write(6,*) ' ---------------------------------'
541     write(6,*) ' Variable: ',cname
542     write(6,*) ' Description: ',cdscrip
543     write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
544     . ' Undefined value = ',undef
545     write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
546     write(6,*) ' Data valid at these times: '
547     ndby3 = ndatebc/3
548     do n = 1, ndby3*3,3
549     write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
550     1000 format(3(2x,i3,':',i8,2x,i8))
551     enddo
552     write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
553 molod 1.23 _END_MASTER( myThid )
554 molod 1.8 endif
555 molod 1.3
556     if( error ) call my_exit (101)
557    
558 molod 1.8 endif
559 molod 1.3
560     C---------- Read SST data if necessary -------------------------------
561    
562 molod 1.5 found = .false.
563     nd = 2
564 molod 1.3
565     c If model time is not within the times of saved sst data
566     c from previous call to getsst then read new data
567    
568 molod 1.23 timemod = dfloat(nymdmod) + dfloat(nhms) /1000000
569     timebc1 = dfloat(nymdbc1) + dfloat(nhmsbc1)/1000000
570     timebc2 = dfloat(nymdbc2) + dfloat(nhmsbc2)/1000000
571    
572 molod 1.5 if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then
573    
574     do while (.not.found .and. nd .le. ndatebc)
575 molod 1.23 timebc2 = dfloat(nymdbc(nd)) + dfloat(nhmsbc(nd))/1000000
576 molod 1.5 if (timebc2 .gt. timemod) then
577     nymdbc1 = nymdbc(nd-1)
578     nymdbc2 = nymdbc(nd)
579     nhmsbc1 = nhmsbc(nd-1)
580     nhmsbc2 = nhmsbc(nd)
581 molod 1.20 call bcdata (iunit,imbc,jmbc,nd,nd+1,sstbc1,sstbc2)
582 molod 1.5 found = .true.
583     else
584     nd = nd + 1
585     endif
586     enddo
587 molod 1.3
588     c Otherwise the data from the last time in getsst surrounds the
589     c current model time.
590    
591 molod 1.5 else
592     found = .true.
593     endif
594 molod 1.3
595 molod 1.5 if (.not.found) then
596     print *, 'STOP: Could not find SST dates for model time.'
597     call my_finalize
598     call my_exit (101)
599     endif
600 molod 1.3
601     C---------- Interpolate SST data ------------------------------------
602    
603 molod 1.5 call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
604     . fac1,fac2)
605 molod 1.3
606 molod 1.5 do j = jm1,jm2
607     do i = im1,im2
608 molod 1.20 sst(i,j,bi,bj) = sstbc1(i+bislot,j+bjslot)*fac1
609     . + sstbc2(i+bislot,j+bjslot)*fac2
610 molod 1.5 enddo
611     enddo
612 molod 1.3
613 molod 1.21
614 molod 1.3 return
615     end
616    
617 molod 1.20 subroutine bcdata (iunit,im,jm,nrec1,nrec2,field1,field2)
618 molod 1.3 C************************************************************************
619     C
620 molod 1.5 C!ROUTINE: BCDATA
621     C!DESCRIPTION: BCDATA reads the data from the file assigned to the
622     C! passed unit number and returns data from the two times
623     C! surrounding the current model time. The two record
624     C! postitions are not assumed to be next to each other.
625 molod 1.3 C
626     C!INPUT PARAMETERS:
627 molod 1.5 C! im number of x points
628     C! im number of x points
629     C! nrec1 record number of the time before the model time
630     C! nrec2 record number of the time after the model time
631 molod 1.3 C
632     C!OUTPUT PARAMETERS:
633 molod 1.20 C! field1(im,jm) data field before the model time
634     C! field2(im,jm) data field after the model time
635 molod 1.3 C
636     C--------------------------------------------------------------------------
637 molod 1.5 implicit none
638 molod 1.3
639 molod 1.20 integer iunit,im,jm,nrec1,nrec2
640 molod 1.3
641 molod 1.20 _RL field1(im,jm)
642     _RL field2(im,jm)
643 molod 1.3
644 molod 1.20 integer i,j
645     real*4 f1(im,jm), f2(im,jm)
646 molod 1.3
647     C--------- Read file -----------------------------------------------
648     read(iunit,rec=nrec1) f1
649     read(iunit,rec=nrec2) f2
650    
651 molod 1.12 #ifdef _BYTESWAPIO
652 molod 1.20 call MDS_BYTESWAPR4( im*jm, f1)
653     call MDS_BYTESWAPR4( im*jm, f2)
654 molod 1.12 #endif
655 molod 1.3 do j=1,jm
656     do i=1,im
657 molod 1.20 field1(i,j) = f1(i,j)
658     field2(i,j) = f2(i,j)
659 molod 1.3 enddo
660     enddo
661    
662     return
663     end
664     subroutine bcheader (iunit, ndmax, nrec,
665 molod 1.20 . cname, cdscrip, im, jm, lat0, lon0, ndatebc,
666 molod 1.5 . nymdbc, nhmsbc, undef, error)
667 molod 1.3 C************************************************************************
668     C
669 molod 1.5 C!ROUTINE: BCHEADER
670     C!DESCRIPTION: BCHEADER reads the header from a file and returns the info.
671 molod 1.3 C
672     C!INPUT PARAMETERS:
673 molod 1.5 C! iunit unit number assigned to the data file
674     C! ndmax maximum number of date/times of the data
675     C! nrec record number of the header info (or assume 1??)
676 molod 1.3 C
677     C!OUTPUT PARAMETERS:
678 molod 1.5 C! cname name of the data in the file header
679     C! cdscrip description of the data in the file header
680     C! im number of x points
681     C! jm number of y points
682     C! lat0 starting latitude for the data grid
683     C! lon0 starting longitude for the data grid
684     C! ndatebc number of date/times of the data in the file
685     C! nymdbc(ndmax) array of dates for the data including century
686     C! nhmsbc(ndmax) array of times for the data
687     C! undef value for undefined values in the data
688     C! error logical TRUE if dataset problem
689 molod 1.3 C
690     C--------------------------------------------------------------------------
691 molod 1.5 implicit none
692 molod 1.3
693 molod 1.5 integer iunit, ndmax, nrec
694 molod 1.3
695 molod 1.5 character*8 cname
696     character*80 cdscrip
697 molod 1.20 character*112 dummy112
698     integer im,jm,ndatebc,nymdbc(ndmax),nhmsbc(ndmax)
699 molod 1.13 _RL lat0,lon0,undef
700 molod 1.5 logical error
701    
702 molod 1.11 integer i
703 molod 1.20 integer*4 im_32,jm_32
704 molod 1.9 integer*4 ndatebc_32,nhmsbc_32(ndmax),nymdbc_32(ndmax)
705 molod 1.5 real*4 lat0_32,lon0_32,undef_32
706 molod 1.3
707     C--------- Read file -----------------------------------------------
708    
709     read(iunit,rec=nrec,err=500) cname, cdscrip,
710 molod 1.20 . im_32, jm_32, lat0_32, lon0_32,
711     . ndatebc_32, undef_32
712    
713 molod 1.12 #ifdef _BYTESWAPIO
714 molod 1.18 call MDS_BYTESWAPI4( 1, im_32)
715     call MDS_BYTESWAPI4( 1, jm_32)
716 molod 1.12 call MDS_BYTESWAPR4( 1, lat0_32)
717     call MDS_BYTESWAPR4( 1, lon0_32)
718 molod 1.20 call MDS_BYTESWAPI4( 1, ndatebc_32)
719 molod 1.12 call MDS_BYTESWAPR4( 1, undef_32)
720     #endif
721 molod 1.5
722 molod 1.20 read(iunit,rec=nrec,err=500) dummy112,
723     . (nymdbc_32(i), nhmsbc_32(i), i=1,ndatebc_32)
724    
725 molod 1.5 im = im_32
726     jm = jm_32
727     lat0 = lat0_32
728     lon0 = lon0_32
729 molod 1.3 undef = undef_32
730    
731     ndatebc = ndatebc_32
732     do i=1,ndatebc
733 molod 1.20 #ifdef _BYTESWAPIO
734     call MDS_BYTESWAPI4( 1, nymdbc_32(i))
735     call MDS_BYTESWAPI4( 1, nhmsbc_32(i))
736     #endif
737 molod 1.3 nymdbc(i) = nymdbc_32(i)
738     nhmsbc(i) = nhmsbc_32(i)
739     enddo
740    
741     return
742     500 continue
743     print *, 'Error reading boundary condition from unit ',iunit
744     error = .true.
745     return
746     end

  ViewVC Help
Powered by ViewVC 1.1.22