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

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

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

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