/[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.25 - (hide annotations) (download)
Mon Mar 12 21:59:10 2007 UTC (17 years, 3 months ago) by molod
Branch: MAIN
Changes since 1.24: +11 -9 lines
Move info around some more

1 molod 1.25 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/update_ocean_exports.F,v 1.24 2007/03/12 21:53:14 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.15 character*40 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.16 sicedata='sice19232.weekly.clim'
198    
199 molod 1.3 c initialize so that first time through they have values for the check
200 molod 1.23 c these values make the iyear .ne. iyearbc true anyways for
201 molod 1.3 c for the first time so first isnt checked below.
202    
203     if (first) then
204     nymdbc(2) = 0
205     nymdbc1 = 0
206     nymdbc2 = 0
207     nhmsbc1 = 0
208     nhmsbc2 = 0
209 molod 1.23 first = .false.
210 molod 1.3 endif
211    
212     C---------- Read in Header file ----------------------------------
213    
214     iyear = nymdmod/10000
215     iyearbc = nymdbc(2)/10000
216    
217     if( iyear.ne.iyearbc ) then
218    
219 molod 1.23 close(iunit)
220 molod 1.15 open (iunit,file=sicedata,form='unformatted',access='direct',
221 molod 1.20 . recl=xsize*ysize*4)
222 molod 1.5 nrec = 1
223     call bcheader (iunit, ndmax, nrec,
224 molod 1.20 . cname, cdscrip, imbc, jmbc, lat0, lon0,
225 molod 1.5 . ndatebc, nymdbc, nhmsbc, undef, error)
226 molod 1.3
227     C--------- Check data for Compatibility ------------------------------
228    
229     C Check for correct data in boundary condition file
230     if (.not.error .and. cname.ne.'SICE') then
231 molod 1.4 write(6,*)'Wrong data in SICE boundary condition file => ',cname
232     error = .true.
233 molod 1.3 endif
234    
235     C Check Horizontal Resolution
236 molod 1.20 if(.not.error.and.imbc*jmbc.ne.xsize*ysize)then
237 molod 1.4 write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
238 molod 1.20 write(6,*) ' B.C. Resolution: ',imbc*jmbc
239     write(6,*) 'Model Resolution: ',xsize*ysize
240 molod 1.4 error = .true.
241 molod 1.3 endif
242    
243     C Check Year
244     iyearbc = nymdbc(2)/10000
245     if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
246 molod 1.5 write(6,*)' B.C. Year DOES NOT match REQUESTED Year!'
247     write(6,*)' B.C. Year: ', iyearbc
248     write(6,*)'Requested Year: ', iyear
249 molod 1.4 error = .true.
250 molod 1.3 endif
251    
252 molod 1.4 if (.not.error) then
253 molod 1.3 C if climatology, fill dates for data with current model year
254 molod 1.4 if (iyearbc.eq.0) then
255     write(6,*)
256     write(6,*) 'Climatological Dataset is being used.'
257     write(6,*) 'Current model year to be used to fill Header Dates'
258     do n = 2, ndatebc-1
259     nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
260     enddo
261 molod 1.3 C For the first date subtract 1 year from the current model NYMD
262 molod 1.4 n = 1
263     nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
264 molod 1.3 C For the last date add 1 year to the current model NYMD
265 molod 1.4 n = ndatebc
266     nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
267     endif
268 molod 1.3
269     C Write out header info
270 molod 1.23 _BEGIN_MASTER( myThid )
271 molod 1.4 write(6,*) ' Updated boundary condition data'
272     write(6,*) ' ---------------------------------'
273     write(6,*) ' Variable: ',cname
274     write(6,*) ' Description: ',cdscrip
275     write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
276     . ' Undefined value = ',undef
277     write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
278     write(6,*) ' Data valid at these times: '
279     ndby3 = ndatebc/3
280     do n = 1, ndby3*3,3
281     write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
282     1000 format(3(2x,i3,':',i8,2x,i8))
283     enddo
284     write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
285 molod 1.23 _END_MASTER( myThid )
286 molod 1.4 endif
287 molod 1.3
288 molod 1.4 endif
289 molod 1.3
290     C---------- Read sice data if necessary -------------------------------
291    
292 molod 1.4 found = .false.
293     nd = 2
294 molod 1.3
295     c If model time is not within the times of saved sice data
296     c from previous call to getsice then read new data
297    
298 molod 1.23 timemod = dfloat(nymdmod) + dfloat(nhms) /1000000
299     timebc1 = dfloat(nymdbc1) + dfloat(nhmsbc1)/1000000
300     timebc2 = dfloat(nymdbc2) + dfloat(nhmsbc2)/1000000
301 molod 1.4
302     if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then
303    
304     do while (.not.found .and. nd .le. ndatebc)
305 molod 1.23 timebc2 = dfloat(nymdbc(nd)) + dfloat(nhmsbc(nd))/1000000
306 molod 1.4 if (timebc2 .gt. timemod) then
307 molod 1.5 nymdbc1 = nymdbc(nd-1)
308     nymdbc2 = nymdbc(nd)
309     nhmsbc1 = nhmsbc(nd-1)
310     nhmsbc2 = nhmsbc(nd)
311 molod 1.20 call bcdata (iunit,imbc,jmbc,nd,nd+1,sicebc1,sicebc2)
312 molod 1.5 found = .true.
313 molod 1.4 else
314 molod 1.5 nd = nd + 1
315 molod 1.4 endif
316     enddo
317 molod 1.3
318     c Otherwise the data from the last time in getsice surrounds the
319     c current model time.
320    
321 molod 1.4 else
322     found = .true.
323     endif
324    
325     if (.not.found) then
326     print *, 'STOP: Could not find SICE dates for model time.'
327     call my_finalize
328     call my_exit (101)
329     endif
330 molod 1.3
331     C---------- Interpolate sice data ------------------------------------
332    
333 molod 1.4 call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
334     . fac1,fac2)
335 molod 1.3
336 molod 1.4 do j = jm1,jm2
337     do i = im1,im2
338 molod 1.20 sice(i,j,bi,bj) = sicebc1(i+bislot,j+bjslot)*fac1
339     . + sicebc2(i+bislot,j+bjslot)*fac2
340 molod 1.3 c average to 0 or 1
341     c -----------------
342 molod 1.4 if (sice(i,j,bi,bj) .ge. 0.5) then
343     sice(i,j,bi,bj) = 1.
344     else
345     sice(i,j,bi,bj) = 0.
346     endif
347     enddo
348     enddo
349 molod 1.3
350     C---------- Fill sice with depth of ice ------------------------------------
351 molod 1.4 do j = jm1,jm2
352     do i = im1,im2
353     if (sice(i,j,bi,bj) .eq. 1.) then
354     sice(i,j,bi,bj) = 3.
355     endif
356     enddo
357     enddo
358 molod 1.3 C---------------------------------------------------------------------------
359    
360     return
361     end
362 molod 1.25 subroutine getsst(iunit,clim,idim1,idim2,jdim1,jdim2,im1,im2,
363     . jm1,jm2,nSumx,nSumy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
364     . sstbc1,sstbc2,first,nymdbc1,nymdbc2,nhmsbc1,nhmsbc2,
365     . nymdbc,nhmsbc,sst,mythid)
366 molod 1.20 C***********************************************************************
367 molod 1.3 C
368     C!ROUTINE: GETSST
369     C!DESCRIPTION: GETSST gets the SST data.
370     C! This routine is adaptable for any frequency
371     C! data upto a daily frequency.
372     C! note: for diurnal data ndmax should be increased.
373     C
374     C!INPUT PARAMETERS:
375 molod 1.5 C! iunit Unit number assigned to the sice data file
376     C! idim1 Start dimension in x-direction
377     C! idim2 End dimension in x-direction
378     C! jdim1 Start dimension in y-direction
379     C! jdim2 End dimension in y-direction
380 molod 1.20 C! im1 Begin of x-direction span for filling sst
381     C! im2 End of x-direction span for filling sst
382     C! jm1 Begin of y-direction span for filling sst
383     C! jm2 End of y-direction span for filling sst
384 molod 1.11 C! nSumx Number of processors in x-direction (local processor)
385     C! nSumy Number of processors in y-direction (local processor)
386 molod 1.20 C! xsize x-dimension of global array
387     C! ysize y-dimension of global array
388 molod 1.5 C! bi Processor number in x-direction (local to processor)
389     C! bj Processor number in y-direction (local to processor)
390 molod 1.20 C! bislot Slot number into global array in x-direction (global)
391     C! bjslot Slot number into global array in y-direction (global)
392 molod 1.5 C! nymd YYMMDD of the current model timestep
393     C! nhms HHMMSS of the model time
394 molod 1.3 C
395     C!OUTPUT PARAMETERS:
396 molod 1.11 C! sst(idim1:idim2,jdim1:jdim2,nSumx,nSumy) Sea surface temperature (K)
397 molod 1.3 C
398     C!ROUTINES CALLED:
399     C
400     C! bcdata Reads the data for a given unit number
401     C! bcheader Reads the header info for a given unit number
402 molod 1.20 C! interp_time Returns weights for linear interpolation
403 molod 1.3 C
404     C--------------------------------------------------------------------------
405    
406     implicit none
407 molod 1.11 #include "SIZE.h"
408 molod 1.5
409 molod 1.11 integer iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSumx,nSumy
410 molod 1.23 integer xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,mythid
411 molod 1.25 logical clim
412 molod 1.3
413 molod 1.20 _RL sstbc1(xsize,ysize)
414     _RL sstbc2(xsize,ysize)
415 molod 1.11 _RL sst(idim1:idim2,jdim1:jdim2,nSumx,nSumy)
416 molod 1.23 integer nhmsbc1,nhmsbc2,nymdbc1,nymdbc2
417     logical first
418 molod 1.3
419 molod 1.5 C Maximum number of dates in one year for the data
420     integer ndmax
421 molod 1.3 parameter (ndmax = 370)
422 molod 1.23 integer nymdbc(ndmax),nhmsbc(ndmax)
423 molod 1.3
424 molod 1.5 character*8 cname
425     character*80 cdscrip
426 molod 1.17 character*20 sstdata
427 molod 1.13 _RL fac1, fac2, lat0, lon0, timebc1, timebc2, timemod, undef
428 molod 1.23 logical found, error
429 molod 1.20 integer i,j,n,nn,iyear,iyearbc,nd,ndby3,imbc,jmbc
430 molod 1.23 integer ndatebc,nrec
431     integer nymdmod
432 molod 1.5
433 molod 1.3
434     C--------- Variable Initialization ---------------------------------
435    
436     data error /.false./
437    
438     c save header info
439 molod 1.23 save imbc,jmbc,lat0,lon0,ndatebc,undef
440 molod 1.3
441     c this only works for between 1950-2050
442     if (nymd .lt. 500101) then
443     nymdmod = 20000000 + nymd
444     else if (nymd .le. 991231) then
445     nymdmod = 19000000 + nymd
446     else
447     nymdmod = nymd
448     endif
449    
450 molod 1.16 sstdata='sst19232.weekly.clim'
451    
452 molod 1.3 c initialize so that first time through they have values for the check
453     c these vaules make the iyear .ne. iyearbc true anyways for
454     c for the first time so first isnt checked below.
455     if (first) then
456     nymdbc(2) = 0
457     nymdbc1 = 0
458     nymdbc2 = 0
459     nhmsbc1 = 0
460     nhmsbc2 = 0
461 molod 1.23 first = .false.
462 molod 1.3 endif
463    
464     C---------- Read in Header file ----------------------------------
465    
466     iyear = nymdmod/10000
467     iyearbc = nymdbc(2)/10000
468    
469     if( iyear.ne.iyearbc ) then
470    
471 molod 1.23 close(iunit)
472 molod 1.15 open (iunit,file=sstdata,form='unformatted',access='direct',
473 molod 1.20 . recl=xsize*ysize*4)
474 molod 1.5 nrec = 1
475     call bcheader (iunit, ndmax, nrec,
476 molod 1.20 . cname, cdscrip, imbc, jmbc, lat0, lon0,
477 molod 1.5 . ndatebc, nymdbc, nhmsbc, undef, error)
478 molod 1.3
479     C--------- Check data for Compatibility
480    
481     C Check for correct data in boundary condition file
482     if (.not.error .and. cname.ne.'SST') then
483 molod 1.5 write(6,*)'Wrong data in SST boundary condition file => ',cname
484     error = .true.
485 molod 1.3 endif
486    
487     C Check Horizontal Resolution
488 molod 1.20 if(.not.error.and.imbc*jmbc.ne.xsize*ysize)then
489 molod 1.5 write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
490 molod 1.20 write(6,*) ' B.C. Resolution: ',imbc*jmbc
491     write(6,*) 'Model Resolution: ',xsize*ysize
492 molod 1.5 error = .true.
493 molod 1.3 endif
494    
495     C Check Year
496     iyearbc = nymdbc(2)/10000
497     if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
498 molod 1.5 write(6,*)' B.C. Year DOES NOT match REQUESTED Year!'
499     write(6,*)' B.C. Year: ', iyearbc
500     write(6,*)'Requested Year: ', iyear
501     error = .true.
502 molod 1.3 endif
503    
504     if (.not.error) then
505     C if climatology, fill dates for data with current model year
506 molod 1.5 if (iyearbc.eq.0) then
507     write(6,*)
508 molod 1.9 write(6,*)'Climatological Dataset is being used.'
509     write(6,*)'Current model year is used to fill Header Dates'
510 molod 1.5 do n = 2, ndatebc-1
511     nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
512     enddo
513 molod 1.3 C For the first date subtract 1 year from the current model NYMD
514 molod 1.5 n = 1
515     nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
516 molod 1.3 C For the last date add 1 year to the current model NYMD
517 molod 1.5 n = ndatebc
518     nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
519     endif
520 molod 1.3
521     C Write out header info
522 molod 1.23 _BEGIN_MASTER( myThid )
523 molod 1.5 write(6,*) ' Updated boundary condition data'
524     write(6,*) ' ---------------------------------'
525     write(6,*) ' Variable: ',cname
526     write(6,*) ' Description: ',cdscrip
527     write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
528     . ' Undefined value = ',undef
529     write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
530     write(6,*) ' Data valid at these times: '
531     ndby3 = ndatebc/3
532     do n = 1, ndby3*3,3
533     write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
534     1000 format(3(2x,i3,':',i8,2x,i8))
535     enddo
536     write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
537 molod 1.23 _END_MASTER( myThid )
538 molod 1.8 endif
539 molod 1.3
540     if( error ) call my_exit (101)
541    
542 molod 1.8 endif
543 molod 1.3
544     C---------- Read SST data if necessary -------------------------------
545    
546 molod 1.5 found = .false.
547     nd = 2
548 molod 1.3
549     c If model time is not within the times of saved sst data
550     c from previous call to getsst then read new data
551    
552 molod 1.23 timemod = dfloat(nymdmod) + dfloat(nhms) /1000000
553     timebc1 = dfloat(nymdbc1) + dfloat(nhmsbc1)/1000000
554     timebc2 = dfloat(nymdbc2) + dfloat(nhmsbc2)/1000000
555    
556 molod 1.5 if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then
557    
558     do while (.not.found .and. nd .le. ndatebc)
559 molod 1.23 timebc2 = dfloat(nymdbc(nd)) + dfloat(nhmsbc(nd))/1000000
560 molod 1.5 if (timebc2 .gt. timemod) then
561     nymdbc1 = nymdbc(nd-1)
562     nymdbc2 = nymdbc(nd)
563     nhmsbc1 = nhmsbc(nd-1)
564     nhmsbc2 = nhmsbc(nd)
565 molod 1.20 call bcdata (iunit,imbc,jmbc,nd,nd+1,sstbc1,sstbc2)
566 molod 1.5 found = .true.
567     else
568     nd = nd + 1
569     endif
570     enddo
571 molod 1.3
572     c Otherwise the data from the last time in getsst surrounds the
573     c current model time.
574    
575 molod 1.5 else
576     found = .true.
577     endif
578 molod 1.3
579 molod 1.5 if (.not.found) then
580     print *, 'STOP: Could not find SST dates for model time.'
581     call my_finalize
582     call my_exit (101)
583     endif
584 molod 1.3
585     C---------- Interpolate SST data ------------------------------------
586    
587 molod 1.5 call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
588     . fac1,fac2)
589 molod 1.3
590 molod 1.5 do j = jm1,jm2
591     do i = im1,im2
592 molod 1.20 sst(i,j,bi,bj) = sstbc1(i+bislot,j+bjslot)*fac1
593     . + sstbc2(i+bislot,j+bjslot)*fac2
594 molod 1.5 enddo
595     enddo
596 molod 1.3
597 molod 1.21
598 molod 1.3 return
599     end
600    
601 molod 1.20 subroutine bcdata (iunit,im,jm,nrec1,nrec2,field1,field2)
602 molod 1.3 C************************************************************************
603     C
604 molod 1.5 C!ROUTINE: BCDATA
605     C!DESCRIPTION: BCDATA reads the data from the file assigned to the
606     C! passed unit number and returns data from the two times
607     C! surrounding the current model time. The two record
608     C! postitions are not assumed to be next to each other.
609 molod 1.3 C
610     C!INPUT PARAMETERS:
611 molod 1.5 C! im number of x points
612     C! im number of x points
613     C! nrec1 record number of the time before the model time
614     C! nrec2 record number of the time after the model time
615 molod 1.3 C
616     C!OUTPUT PARAMETERS:
617 molod 1.20 C! field1(im,jm) data field before the model time
618     C! field2(im,jm) data field after the model time
619 molod 1.3 C
620     C--------------------------------------------------------------------------
621 molod 1.5 implicit none
622 molod 1.3
623 molod 1.20 integer iunit,im,jm,nrec1,nrec2
624 molod 1.3
625 molod 1.20 _RL field1(im,jm)
626     _RL field2(im,jm)
627 molod 1.3
628 molod 1.20 integer i,j
629     real*4 f1(im,jm), f2(im,jm)
630 molod 1.3
631     C--------- Read file -----------------------------------------------
632     read(iunit,rec=nrec1) f1
633     read(iunit,rec=nrec2) f2
634    
635 molod 1.12 #ifdef _BYTESWAPIO
636 molod 1.20 call MDS_BYTESWAPR4( im*jm, f1)
637     call MDS_BYTESWAPR4( im*jm, f2)
638 molod 1.12 #endif
639 molod 1.3 do j=1,jm
640     do i=1,im
641 molod 1.20 field1(i,j) = f1(i,j)
642     field2(i,j) = f2(i,j)
643 molod 1.3 enddo
644     enddo
645    
646     return
647     end
648     subroutine bcheader (iunit, ndmax, nrec,
649 molod 1.20 . cname, cdscrip, im, jm, lat0, lon0, ndatebc,
650 molod 1.5 . nymdbc, nhmsbc, undef, error)
651 molod 1.3 C************************************************************************
652     C
653 molod 1.5 C!ROUTINE: BCHEADER
654     C!DESCRIPTION: BCHEADER reads the header from a file and returns the info.
655 molod 1.3 C
656     C!INPUT PARAMETERS:
657 molod 1.5 C! iunit unit number assigned to the data file
658     C! ndmax maximum number of date/times of the data
659     C! nrec record number of the header info (or assume 1??)
660 molod 1.3 C
661     C!OUTPUT PARAMETERS:
662 molod 1.5 C! cname name of the data in the file header
663     C! cdscrip description of the data in the file header
664     C! im number of x points
665     C! jm number of y points
666     C! lat0 starting latitude for the data grid
667     C! lon0 starting longitude for the data grid
668     C! ndatebc number of date/times of the data in the file
669     C! nymdbc(ndmax) array of dates for the data including century
670     C! nhmsbc(ndmax) array of times for the data
671     C! undef value for undefined values in the data
672     C! error logical TRUE if dataset problem
673 molod 1.3 C
674     C--------------------------------------------------------------------------
675 molod 1.5 implicit none
676 molod 1.3
677 molod 1.5 integer iunit, ndmax, nrec
678 molod 1.3
679 molod 1.5 character*8 cname
680     character*80 cdscrip
681 molod 1.20 character*112 dummy112
682     integer im,jm,ndatebc,nymdbc(ndmax),nhmsbc(ndmax)
683 molod 1.13 _RL lat0,lon0,undef
684 molod 1.5 logical error
685    
686 molod 1.11 integer i
687 molod 1.20 integer*4 im_32,jm_32
688 molod 1.9 integer*4 ndatebc_32,nhmsbc_32(ndmax),nymdbc_32(ndmax)
689 molod 1.5 real*4 lat0_32,lon0_32,undef_32
690 molod 1.3
691     C--------- Read file -----------------------------------------------
692    
693     read(iunit,rec=nrec,err=500) cname, cdscrip,
694 molod 1.20 . im_32, jm_32, lat0_32, lon0_32,
695     . ndatebc_32, undef_32
696    
697 molod 1.12 #ifdef _BYTESWAPIO
698 molod 1.18 call MDS_BYTESWAPI4( 1, im_32)
699     call MDS_BYTESWAPI4( 1, jm_32)
700 molod 1.12 call MDS_BYTESWAPR4( 1, lat0_32)
701     call MDS_BYTESWAPR4( 1, lon0_32)
702 molod 1.20 call MDS_BYTESWAPI4( 1, ndatebc_32)
703 molod 1.12 call MDS_BYTESWAPR4( 1, undef_32)
704     #endif
705 molod 1.5
706 molod 1.20 read(iunit,rec=nrec,err=500) dummy112,
707     . (nymdbc_32(i), nhmsbc_32(i), i=1,ndatebc_32)
708    
709 molod 1.5 im = im_32
710     jm = jm_32
711     lat0 = lat0_32
712     lon0 = lon0_32
713 molod 1.3 undef = undef_32
714    
715     ndatebc = ndatebc_32
716     do i=1,ndatebc
717 molod 1.20 #ifdef _BYTESWAPIO
718     call MDS_BYTESWAPI4( 1, nymdbc_32(i))
719     call MDS_BYTESWAPI4( 1, nhmsbc_32(i))
720     #endif
721 molod 1.3 nymdbc(i) = nymdbc_32(i)
722     nhmsbc(i) = nhmsbc_32(i)
723     enddo
724    
725     return
726     500 continue
727     print *, 'Error reading boundary condition from unit ',iunit
728     error = .true.
729     return
730     end

  ViewVC Help
Powered by ViewVC 1.1.22