/[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.33 - (hide annotations) (download)
Wed Sep 22 22:21:41 2010 UTC (13 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint65, checkpoint63, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, HEAD
Changes since 1.32: +2 -2 lines
remove "??)" in comments (trigger a warning from xmakedepend)

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

  ViewVC Help
Powered by ViewVC 1.1.22