Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:40:45 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
e337e4ca8c Andr*0002        subroutine update_ocean_exports (myTime, myIter, myThid)
                0003 c----------------------------------------------------------------------
                0004 c  Subroutine update_ocean_exports - 'Wrapper' routine to update
3daafce20b Jean*0005 c        the fields related to the ocean surface that are needed
e17ffcecb1 Andr*0006 c        by fizhi (sst and sea ice extent).
e337e4ca8c Andr*0007 c
                0008 c Call:  getsst  (Return the current sst field-read dataset if needed)
                0009 c        getsice (Return the current sea ice field-read data if needed)
                0010 c-----------------------------------------------------------------------
                0011        implicit none
                0012 #include "SIZE.h"
                0013 #include "GRID.h"
f4a0368053 Andr*0014 #include "fizhi_ocean_coms.h"
e337e4ca8c Andr*0015 #include "EEPARAMS.h"
e17ffcecb1 Andr*0016 #include "chronos.h"
819a3c957e Andr*0017 #ifdef ALLOW_EXCH2
f9f661930b Jean*0018 #include "W2_EXCH2_SIZE.h"
819a3c957e Andr*0019 #include "W2_EXCH2_TOPOLOGY.h"
                0020 #endif /* ALLOW_EXCH2 */
e337e4ca8c Andr*0021 
3768927683 Andr*0022        integer myIter, myThid
                0023        _RL myTime
e337e4ca8c Andr*0024 
a6b4f97600 Jean*0025        INTEGER xySize
                0026 #if defined(ALLOW_EXCH2)
                0027        PARAMETER ( xySize = W2_ioBufferSize )
                0028 #else
                0029        PARAMETER ( xySize = Nx*Ny )
                0030 #endif
819a3c957e Andr*0031        integer i, j, bi, bj, bislot, bjslot
e337e4ca8c Andr*0032        integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
819a3c957e Andr*0033        integer xsize, ysize
a6b4f97600 Jean*0034        _RL        sstmin
e17ffcecb1 Andr*0035        parameter ( sstmin = 273.16 )
e337e4ca8c Andr*0036 
a6b4f97600 Jean*0037        _RL sst1 (xySize), sst2 (xySize)
                0038        _RL sice1(xySize), sice2(xySize)
                0039 c      _RL sst1(xsize,ysize),sst2(xsize,ysize)
                0040 c      _RL sice1(xsize,ysize),sice2(xsize,ysize)
3f946231fb Andr*0041        integer nymd1sst(nSx,nSy),nymd2sst(nSx,nSy)
                0042        integer nymd1sice(nSx,nSy),nymd2sice(nSx,nSy)
                0043        integer nhms1sst(nSx,nSy),nhms2sst(nSx,nSy)
                0044        integer nhms1sice(nSx,nSy),nhms2sice(nSx,nSy)
                0045        integer sstdates(370,nSx,nSy),sicedates(370,nSx,nSy)
                0046        integer ssttimes(370,nSx,nSy),sicetimes(370,nSx,nSy)
                0047        logical first(nSx,nSy)
                0048        integer nSxnSy
                0049        parameter(nSxnSy = nSx*nSy)
                0050        data first/nSxnSy*.true./
                0051 
                0052        save nymd1sst,nymd2sst,nymd1sice,nymd2sice
                0053        save nhms1sst,nhms2sst,nhms1sice,nhms2sice
819a3c957e Andr*0054        save sst1, sst2, sice1, sice2
3f946231fb Andr*0055        save sstdates, sicedates
                0056        save ssttimes, sicetimes
819a3c957e Andr*0057 
a6b4f97600 Jean*0058 #if defined(ALLOW_EXCH2)
                0059        xsize = exch2_global_Nx
                0060        ysize = exch2_global_Ny
                0061 #else
                0062        xsize = Nx
                0063        ysize = Ny
                0064 #endif
218575850c Andr*0065        idim1 = 1-OLx
                0066        idim2 = sNx+OLx
                0067        jdim1 = 1-OLy
                0068        jdim2 = sNy+OLy
                0069        im1 = 1
                0070        im2 = sNx
                0071        jm1 = 1
                0072        jm2 = sNy
                0073 
                0074 C***********************************************************************
                0075 
                0076        DO BJ = myByLo(myThid),myByHi(myThid)
                0077        DO BI = myBxLo(myThid),myBxHi(myThid)
819a3c957e Andr*0078 #if defined(ALLOW_EXCH2)
c424ee7cc7 Jean*0079        bislot = exch2_txglobalo(W2_myTileList(bi,bj))-1
                0080        bjslot = exch2_tyglobalo(W2_myTileList(bi,bj))-1
819a3c957e Andr*0081 #else
                0082        bislot = myXGlobalLo-1+(bi-1)*sNx
                0083        bjslot = myYGlobalLo-1+(bj-1)*sNy
                0084 #endif
e17ffcecb1 Andr*0085 
8a8ee991ea Andr*0086        call getsst(ksst,sstclim,idim1,idim2,jdim1,jdim2,im1,im2,
                0087      .  jm1,jm2,nSx,nSy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
                0088      .  sst1,sst2,first(bi,bj),nymd1sst(bi,bj),nymd2sst(bi,bj),
3f946231fb Andr*0089      .  nhms1sst(bi,bj),nhms2sst(bi,bj),sstdates(1,bi,bj),
                0090      .  ssttimes(1,bi,bj),sst,myThid)
8a8ee991ea Andr*0091        call getsice(kice,siceclim,idim1,idim2,jdim1,jdim2,im1,im2,
                0092      .  jm1,jm2,nSx,nSy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
                0093      .  sice1,sice2,first(bi,bj),nymd1sice(bi,bj),nymd2sice(bi,bj),
3f946231fb Andr*0094      .  nhms1sice(bi,bj),nhms2sice(bi,bj),sicedates(1,bi,bj),
                0095      .  sicetimes(1,bi,bj),sice,myThid)
e17ffcecb1 Andr*0096 
                0097 c Check for Minimum Open-Water SST
                0098 c --------------------------------
218575850c Andr*0099        do j=jm1,jm2
                0100        do i=im1,im2
                0101        if(sice(i,j,bi,bj).eq.0.0 .and. sst(i,j,bi,bj).lt.sstmin)
                0102      .                                          sst(i,j,bi,bj) = sstmin
e337e4ca8c Andr*0103        enddo
                0104        enddo
                0105 
218575850c Andr*0106        ENDDO
                0107        ENDDO
6637358eea Jean*0108        _EXCH_XY_RL(sst,myThid)
                0109        _EXCH_XY_RL(sice,myThid)
218575850c Andr*0110 
e337e4ca8c Andr*0111        return
                0112        end
e17ffcecb1 Andr*0113 
8d4247b0d2 Andr*0114        subroutine getsice(iunit,clim,idim1,idim2,jdim1,jdim2,im1,im2,
                0115      .   jm1,jm2,nSumx,nSumy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
                0116      .   sicebc1,sicebc2,first,nymdbc1,nymdbc2,nhmsbc1,nhmsbc2,
                0117      .   nymdbc,nhmsbc,sice,mythid)
819a3c957e Andr*0118 C***********************************************************************
e17ffcecb1 Andr*0119 C
a6b4f97600 Jean*0120 C!ROUTINE: GETSICE
                0121 C!DESCRIPTION:  GETSICE returns the sea ice depth.
                0122 C!              This routine is adaptable for any frequency
                0123 C!              data upto a daily frequency.
                0124 C!              note: for diurnal data ndmax should be increased.
e17ffcecb1 Andr*0125 C
                0126 C!INPUT PARAMETERS:
a6b4f97600 Jean*0127 C!      iunit     Unit number assigned to the sice data file
218575850c Andr*0128 C!      idim1     Start dimension in x-direction
                0129 C!      idim2     End dimension in x-direction
                0130 C!      jdim1     Start dimension in y-direction
                0131 C!      jdim2     End dimension in y-direction
                0132 C!      im1       Begin of x-direction span for filling sice
                0133 C!      im2       End of x-direction span for filling sice
                0134 C!      jm1       Begin of y-direction span for filling sice
                0135 C!      jm2       End of y-direction span for filling sice
4f001d3c64 Andr*0136 C!      nSumx     Number of processors in x-direction (local processor)
                0137 C!      nSumy     Number of processors in y-direction (local processor)
819a3c957e Andr*0138 C!      xsize      Number of processors in x-direction (global)
                0139 C!      ysize      Number of processors in y-direction (global)
218575850c Andr*0140 C!      bi        Processor number in x-direction (local to processor)
                0141 C!      bj        Processor number in y-direction (local to processor)
819a3c957e Andr*0142 C!      bislot  Processor number in x-direction (global)
                0143 C!      bjslot  Processor number in y-direction (global)
a6b4f97600 Jean*0144 C!      nymd    YYMMDD of the current model timestep
                0145 C!      nhms    HHMMSS of the model time
e17ffcecb1 Andr*0146 C
                0147 C!OUTPUT PARAMETERS:
a6b4f97600 Jean*0148 C!      sice(idim1:idim2,jdim1:jdim2,nSumx,nSumy) Sea ice depth in meters
e17ffcecb1 Andr*0149 C
                0150 C!ROUTINES CALLED:
                0151 C
a6b4f97600 Jean*0152 C!      bcdata       Reads the data for a given unit number
                0153 C!      bcheader     Reads the header info for a given unit number
218575850c Andr*0154 C!      interp_time  Returns weights for linear interpolation
e17ffcecb1 Andr*0155 C
                0156 C--------------------------------------------------------------------------
                0157 
                0158       implicit none
4f001d3c64 Andr*0159 #include "SIZE.h"
218575850c Andr*0160 
4f001d3c64 Andr*0161       integer iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSumx,nSumy
3f946231fb Andr*0162       integer xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,mythid
8d4247b0d2 Andr*0163       logical clim
218575850c Andr*0164 
819a3c957e Andr*0165       _RL sicebc1(xsize,ysize)
                0166       _RL sicebc2(xsize,ysize)
4f001d3c64 Andr*0167       _RL sice(idim1:idim2,jdim1:jdim2,nSumx,nSumy)
3f946231fb Andr*0168       integer nhmsbc1,nhmsbc2,nymdbc1,nymdbc2
                0169       logical first
e17ffcecb1 Andr*0170 
218575850c Andr*0171 C Maximum number of dates in one year for the data
a6b4f97600 Jean*0172       integer ndmax
e17ffcecb1 Andr*0173       parameter (ndmax = 370)
3f946231fb Andr*0174       integer nymdbc(ndmax),nhmsbc(ndmax)
e17ffcecb1 Andr*0175 
218575850c Andr*0176       character*8  cname
                0177       character*80 cdscrip
e4382a24c8 Andr*0178       character*22 sicedata
a456aa407c Andr*0179       _RL fac1, fac2, lat0, lon0, timebc1, timebc2, timemod, undef
3f946231fb Andr*0180       logical found, error
819a3c957e Andr*0181       integer i,j,n,nn,iyear,iyearbc,nd,ndby3,imbc,jmbc
3f946231fb Andr*0182       integer ndatebc,nrec
                0183       integer nymdmod
e17ffcecb1 Andr*0184 
                0185 
                0186 C--------- Variable Initialization ---------------------------------
                0187 
                0188       data error /.false./
                0189 
                0190 c  save header info
3f946231fb Andr*0191       save imbc,jmbc,lat0,lon0,ndatebc,undef
e17ffcecb1 Andr*0192 
                0193 c  this only works for between 1950-2050
                0194       if (nymd .lt. 500101) then
                0195         nymdmod = 20000000 + nymd
                0196       else if (nymd .le. 991231) then
                0197         nymdmod = 19000000 + nymd
                0198       else
                0199         nymdmod = nymd
                0200       endif
                0201 
e4382a24c8 Andr*0202       iyear   = nymdmod/10000
                0203       if(clim) then
f1edb8ebdb Andr*0204        if(xsize.eq.192)sicedata='sice19232.weekly.clim'
                0205        if(xsize.eq.612)sicedata='sice612102.weekly.clim'
e4382a24c8 Andr*0206       else
f1edb8ebdb Andr*0207        if(xsize.eq.192)
                0208      .           WRITE(sicedata,'(A,I4)')'sice19232.weekly.y',iyear
                0209        if(xsize.eq.612)
                0210      .           WRITE(sicedata,'(A,I4)')'sice612102.weekly.y',iyear
e4382a24c8 Andr*0211       endif
30b199c293 Andr*0212 
e17ffcecb1 Andr*0213 c  initialize so that first time through they have values for the check
3f946231fb Andr*0214 c  these values make the iyear .ne. iyearbc true anyways for
e17ffcecb1 Andr*0215 c  for the first time so first isnt checked below.
                0216 
                0217       if (first) then
                0218         nymdbc(2) = 0
                0219         nymdbc1   = 0
                0220         nymdbc2   = 0
                0221         nhmsbc1   = 0
                0222         nhmsbc2   = 0
3f946231fb Andr*0223         first = .false.
e17ffcecb1 Andr*0224       endif
                0225 
                0226 C---------- Read in Header file ----------------------------------
                0227 
                0228       iyearbc = nymdbc(2)/10000
                0229 
                0230       if( iyear.ne.iyearbc ) then
                0231 
3f946231fb Andr*0232        close(iunit)
8c99ca8c64 Andr*0233        open (iunit,file=sicedata,form='unformatted',access='direct',
819a3c957e Andr*0234      .                                         recl=xsize*ysize*4)
598e83bd78 Andr*0235        nrec = 1
                0236        call bcheader (iunit, ndmax, nrec,
a6b4f97600 Jean*0237      .          cname, cdscrip, imbc, jmbc, lat0, lon0,
598e83bd78 Andr*0238      .          ndatebc, nymdbc, nhmsbc, undef, error)
e17ffcecb1 Andr*0239 
                0240 C--------- Check data for Compatibility ------------------------------
                0241 
                0242 C Check for correct data in boundary condition file
                0243        if (.not.error .and. cname.ne.'SICE') then
218575850c Andr*0244         write(6,*)'Wrong data in SICE boundary condition file => ',cname
                0245         error = .true.
e17ffcecb1 Andr*0246        endif
                0247 
                0248 C Check Horizontal Resolution
819a3c957e Andr*0249        if(.not.error.and.imbc*jmbc.ne.xsize*ysize)then
218575850c Andr*0250         write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
819a3c957e Andr*0251         write(6,*) ' B.C. Resolution:  ',imbc*jmbc
                0252         write(6,*) 'Model Resolution:  ',xsize*ysize
218575850c Andr*0253         error = .true.
e17ffcecb1 Andr*0254        endif
                0255 
a6b4f97600 Jean*0256 C Check Year
e17ffcecb1 Andr*0257        iyearbc = nymdbc(2)/10000
                0258        if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
598e83bd78 Andr*0259         write(6,*)'     B.C. Year DOES NOT match REQUESTED Year!'
                0260         write(6,*)'     B.C. Year:  ', iyearbc
                0261         write(6,*)'Requested Year:  ', iyear
218575850c Andr*0262         error = .true.
e17ffcecb1 Andr*0263        endif
                0264 
218575850c Andr*0265        if (.not.error)   then
e17ffcecb1 Andr*0266 C if climatology, fill dates for data with current model year
a6b4f97600 Jean*0267         if (iyearbc.eq.0) then
218575850c Andr*0268          write(6,*)
a6b4f97600 Jean*0269          write(6,*) 'Climatological Dataset is being used.'
218575850c Andr*0270          write(6,*) 'Current model year to be used to fill Header Dates'
                0271          do n = 2, ndatebc-1
a6b4f97600 Jean*0272           nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
218575850c Andr*0273          enddo
e17ffcecb1 Andr*0274 C  For the first date subtract 1 year from the current model NYMD
218575850c Andr*0275          n = 1
a6b4f97600 Jean*0276          nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
e17ffcecb1 Andr*0277 C  For the last date add 1 year to the current model NYMD
218575850c Andr*0278          n = ndatebc
a6b4f97600 Jean*0279          nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
218575850c Andr*0280         endif
e17ffcecb1 Andr*0281 
                0282 C  Write out header info
3f946231fb Andr*0283         _BEGIN_MASTER( myThid )
218575850c Andr*0284         write(6,*) ' Updated boundary condition data'
                0285         write(6,*) ' ---------------------------------'
                0286         write(6,*) ' Variable: ',cname
                0287         write(6,*) ' Description: ',cdscrip
                0288         write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
                0289      .                                       ' Undefined value = ',undef
                0290         write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
                0291         write(6,*) ' Data valid at these times: '
                0292         ndby3 = ndatebc/3
                0293         do n = 1, ndby3*3,3
                0294          write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
                0295  1000    format(3(2x,i3,':',i8,2x,i8))
                0296         enddo
                0297         write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
3f946231fb Andr*0298         _END_MASTER( myThid )
a6b4f97600 Jean*0299        endif
e17ffcecb1 Andr*0300 
218575850c Andr*0301       endif
e17ffcecb1 Andr*0302 
                0303 C---------- Read sice data if necessary -------------------------------
                0304 
218575850c Andr*0305       found = .false.
                0306       nd = 2
e17ffcecb1 Andr*0307 
a6b4f97600 Jean*0308 c  If model time is not within the times of saved sice data
                0309 c  from previous call to getsice then read new data
e17ffcecb1 Andr*0310 
3f946231fb Andr*0311       timemod = dfloat(nymdmod) + dfloat(nhms)   /1000000
                0312       timebc1 = dfloat(nymdbc1) + dfloat(nhmsbc1)/1000000
                0313       timebc2 = dfloat(nymdbc2) + dfloat(nhmsbc2)/1000000
e17ffcecb1 Andr*0314 
218575850c Andr*0315       if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then
e17ffcecb1 Andr*0316 
218575850c Andr*0317        do while (.not.found .and. nd .le. ndatebc)
3f946231fb Andr*0318         timebc2 = dfloat(nymdbc(nd)) + dfloat(nhmsbc(nd))/1000000
a6b4f97600 Jean*0319         if (timebc2 .gt. timemod) then
598e83bd78 Andr*0320          nymdbc1 = nymdbc(nd-1)
                0321          nymdbc2 = nymdbc(nd)
                0322          nhmsbc1 = nhmsbc(nd-1)
                0323          nhmsbc2 = nhmsbc(nd)
819a3c957e Andr*0324          call bcdata (iunit,imbc,jmbc,nd,nd+1,sicebc1,sicebc2)
598e83bd78 Andr*0325          found = .true.
218575850c Andr*0326         else
598e83bd78 Andr*0327          nd = nd + 1
e17ffcecb1 Andr*0328         endif
218575850c Andr*0329        enddo
                0330 
a6b4f97600 Jean*0331 c  Otherwise the data from the last time in getsice surrounds the
218575850c Andr*0332 c  current model time.
                0333 
                0334       else
                0335        found = .true.
                0336       endif
                0337 
                0338       if (.not.found) then
a6b4f97600 Jean*0339        print *, 'STOP: Could not find SICE dates for model time.'
218575850c Andr*0340        call my_finalize
                0341        call my_exit (101)
                0342       endif
e17ffcecb1 Andr*0343 
                0344 C---------- Interpolate sice data ------------------------------------
                0345 
218575850c Andr*0346       call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
                0347      .                                                       fac1,fac2)
e17ffcecb1 Andr*0348 
218575850c Andr*0349       do j = jm1,jm2
                0350       do i = im1,im2
a6b4f97600 Jean*0351        sice(i,j,bi,bj) = sicebc1(i+bislot,j+bjslot)*fac1
                0352      .                 + sicebc2(i+bislot,j+bjslot)*fac2
e17ffcecb1 Andr*0353 c average to 0 or 1
                0354 c -----------------
218575850c Andr*0355        if (sice(i,j,bi,bj) .ge. 0.5) then
                0356         sice(i,j,bi,bj) = 1.
                0357        else
                0358         sice(i,j,bi,bj) = 0.
                0359        endif
                0360       enddo
                0361       enddo
e17ffcecb1 Andr*0362 
                0363 C---------- Fill sice with depth of ice ------------------------------------
218575850c Andr*0364       do j = jm1,jm2
                0365       do i = im1,im2
a6b4f97600 Jean*0366        if (sice(i,j,bi,bj) .eq. 1.) then
218575850c Andr*0367         sice(i,j,bi,bj) = 3.
                0368        endif
                0369       enddo
                0370       enddo
e17ffcecb1 Andr*0371 C---------------------------------------------------------------------------
a6b4f97600 Jean*0372 
e17ffcecb1 Andr*0373       return
                0374       end
8d4247b0d2 Andr*0375       subroutine getsst(iunit,clim,idim1,idim2,jdim1,jdim2,im1,im2,
                0376      .   jm1,jm2,nSumx,nSumy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
                0377      .   sstbc1,sstbc2,first,nymdbc1,nymdbc2,nhmsbc1,nhmsbc2,
                0378      .   nymdbc,nhmsbc,sst,mythid)
819a3c957e Andr*0379 C***********************************************************************
e17ffcecb1 Andr*0380 C
a6b4f97600 Jean*0381 C!ROUTINE: GETSST
                0382 C!DESCRIPTION:  GETSST gets the SST data.
                0383 C!              This routine is adaptable for any frequency
                0384 C!              data upto a daily frequency.
                0385 C!              note: for diurnal data ndmax should be increased.
e17ffcecb1 Andr*0386 C
                0387 C!INPUT PARAMETERS:
a6b4f97600 Jean*0388 C!      iunit     Unit number assigned to the sice data file
598e83bd78 Andr*0389 C!      idim1     Start dimension in x-direction
                0390 C!      idim2     End dimension in x-direction
                0391 C!      jdim1     Start dimension in y-direction
                0392 C!      jdim2     End dimension in y-direction
819a3c957e Andr*0393 C!      im1       Begin of x-direction span for filling sst
                0394 C!      im2       End of x-direction span for filling sst
                0395 C!      jm1       Begin of y-direction span for filling sst
                0396 C!      jm2       End of y-direction span for filling sst
4f001d3c64 Andr*0397 C!      nSumx     Number of processors in x-direction (local processor)
                0398 C!      nSumy     Number of processors in y-direction (local processor)
819a3c957e Andr*0399 C!      xsize     x-dimension of global array
                0400 C!      ysize     y-dimension of global array
598e83bd78 Andr*0401 C!      bi        Processor number in x-direction (local to processor)
                0402 C!      bj        Processor number in y-direction (local to processor)
819a3c957e Andr*0403 C!      bislot    Slot number into global array in x-direction (global)
                0404 C!      bjslot    Slot number into global array in y-direction (global)
a6b4f97600 Jean*0405 C!      nymd      YYMMDD of the current model timestep
                0406 C!      nhms      HHMMSS of the model time
e17ffcecb1 Andr*0407 C
                0408 C!OUTPUT PARAMETERS:
a6b4f97600 Jean*0409 C!     sst(idim1:idim2,jdim1:jdim2,nSumx,nSumy) Sea surface temperature (K)
e17ffcecb1 Andr*0410 C
                0411 C!ROUTINES CALLED:
                0412 C
a6b4f97600 Jean*0413 C!     bcdata     Reads the data for a given unit number
                0414 C!     bcheader   Reads the header info for a given unit number
                0415 C!     interp_time  Returns weights for linear interpolation
e17ffcecb1 Andr*0416 C
                0417 C--------------------------------------------------------------------------
                0418 
                0419       implicit none
4f001d3c64 Andr*0420 #include "SIZE.h"
e17ffcecb1 Andr*0421 
4f001d3c64 Andr*0422       integer iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSumx,nSumy
3f946231fb Andr*0423       integer xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,mythid
8d4247b0d2 Andr*0424       logical clim
e17ffcecb1 Andr*0425 
819a3c957e Andr*0426       _RL sstbc1(xsize,ysize)
                0427       _RL sstbc2(xsize,ysize)
4f001d3c64 Andr*0428       _RL sst(idim1:idim2,jdim1:jdim2,nSumx,nSumy)
3f946231fb Andr*0429       integer nhmsbc1,nhmsbc2,nymdbc1,nymdbc2
                0430       logical first
598e83bd78 Andr*0431 
                0432 C Maximum number of dates in one year for the data
a6b4f97600 Jean*0433       integer ndmax
e17ffcecb1 Andr*0434       parameter (ndmax = 370)
3f946231fb Andr*0435       integer nymdbc(ndmax),nhmsbc(ndmax)
e17ffcecb1 Andr*0436 
598e83bd78 Andr*0437       character*8  cname
                0438       character*80 cdscrip
e4382a24c8 Andr*0439       character*21 sstdata
a456aa407c Andr*0440       _RL fac1, fac2, lat0, lon0, timebc1, timebc2, timemod, undef
3f946231fb Andr*0441       logical found, error
819a3c957e Andr*0442       integer i,j,n,nn,iyear,iyearbc,nd,ndby3,imbc,jmbc
3f946231fb Andr*0443       integer ndatebc,nrec
                0444       integer nymdmod
598e83bd78 Andr*0445 
e17ffcecb1 Andr*0446 
                0447 C--------- Variable Initialization ---------------------------------
                0448 
                0449       data error /.false./
                0450 
                0451 c  save header info
3f946231fb Andr*0452       save imbc,jmbc,lat0,lon0,ndatebc,undef
e17ffcecb1 Andr*0453 
                0454 c  this only works for between 1950-2050
                0455       if (nymd .lt. 500101) then
                0456         nymdmod = 20000000 + nymd
                0457       else if (nymd .le. 991231) then
                0458         nymdmod = 19000000 + nymd
                0459       else
                0460         nymdmod = nymd
                0461       endif
                0462 
e4382a24c8 Andr*0463       iyear   = nymdmod/10000
                0464       if(clim) then
f1edb8ebdb Andr*0465        if(xsize.eq.192)sstdata='sst19232.weekly.clim'
                0466        if(xsize.eq.612)sstdata='sst612102.weekly.clim'
e4382a24c8 Andr*0467       else
f1edb8ebdb Andr*0468        if(xsize.eq.192)
                0469      .           WRITE(sstdata,'(A,I4)')'sst19232.weekly.y',iyear
                0470        if(xsize.eq.612)
                0471      .           WRITE(sstdata,'(A,I4)')'sst612102.weekly.y',iyear
e4382a24c8 Andr*0472       endif
30b199c293 Andr*0473 
e17ffcecb1 Andr*0474 c  initialize so that first time through they have values for the check
                0475 c  these vaules make the iyear .ne. iyearbc true anyways for
                0476 c  for the first time so first isnt checked below.
                0477       if (first) then
                0478         nymdbc(2) = 0
                0479         nymdbc1   = 0
                0480         nymdbc2   = 0
                0481         nhmsbc1   = 0
                0482         nhmsbc2   = 0
3f946231fb Andr*0483         first = .false.
e17ffcecb1 Andr*0484       endif
                0485 
                0486 C---------- Read in Header file ----------------------------------
                0487 
                0488       iyearbc = nymdbc(2)/10000
                0489 
                0490       if( iyear.ne.iyearbc ) then
                0491 
3f946231fb Andr*0492        close(iunit)
8c99ca8c64 Andr*0493        open (iunit,file=sstdata,form='unformatted',access='direct',
819a3c957e Andr*0494      .                                        recl=xsize*ysize*4)
598e83bd78 Andr*0495        nrec = 1
                0496        call bcheader (iunit, ndmax, nrec,
a6b4f97600 Jean*0497      .          cname, cdscrip, imbc, jmbc, lat0, lon0,
598e83bd78 Andr*0498      .          ndatebc, nymdbc, nhmsbc, undef, error)
e17ffcecb1 Andr*0499 
                0500 C--------- Check data for Compatibility
                0501 
                0502 C Check for correct data in boundary condition file
                0503        if (.not.error .and. cname.ne.'SST') then
598e83bd78 Andr*0504         write(6,*)'Wrong data in SST boundary condition file => ',cname
                0505         error = .true.
e17ffcecb1 Andr*0506        endif
                0507 
                0508 C Check Horizontal Resolution
819a3c957e Andr*0509        if(.not.error.and.imbc*jmbc.ne.xsize*ysize)then
598e83bd78 Andr*0510         write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
819a3c957e Andr*0511         write(6,*) ' B.C. Resolution:  ',imbc*jmbc
                0512         write(6,*) 'Model Resolution:  ',xsize*ysize
598e83bd78 Andr*0513         error = .true.
e17ffcecb1 Andr*0514        endif
                0515 
a6b4f97600 Jean*0516 C Check Year
e17ffcecb1 Andr*0517        iyearbc = nymdbc(2)/10000
                0518        if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
598e83bd78 Andr*0519         write(6,*)'     B.C. Year DOES NOT match REQUESTED Year!'
                0520         write(6,*)'     B.C. Year:  ', iyearbc
                0521         write(6,*)'Requested Year:  ', iyear
                0522         error = .true.
e17ffcecb1 Andr*0523        endif
                0524 
                0525        if (.not.error)   then
                0526 C if climatology, fill dates for data with current model year
a6b4f97600 Jean*0527         if (iyearbc.eq.0) then
598e83bd78 Andr*0528          write(6,*)
a6b4f97600 Jean*0529          write(6,*)'Climatological Dataset is being used.'
613fa3996d Andr*0530          write(6,*)'Current model year is used to fill Header Dates'
598e83bd78 Andr*0531          do n = 2, ndatebc-1
a6b4f97600 Jean*0532           nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
598e83bd78 Andr*0533          enddo
e17ffcecb1 Andr*0534 C  For the first date subtract 1 year from the current model NYMD
598e83bd78 Andr*0535          n = 1
a6b4f97600 Jean*0536          nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
e17ffcecb1 Andr*0537 C  For the last date add 1 year to the current model NYMD
598e83bd78 Andr*0538          n = ndatebc
a6b4f97600 Jean*0539          nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
598e83bd78 Andr*0540         endif
e17ffcecb1 Andr*0541 
                0542 C  Write out header info
3f946231fb Andr*0543         _BEGIN_MASTER( myThid )
598e83bd78 Andr*0544         write(6,*) ' Updated boundary condition data'
                0545         write(6,*) ' ---------------------------------'
                0546         write(6,*) ' Variable: ',cname
                0547         write(6,*) ' Description: ',cdscrip
                0548         write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
                0549      .                                       ' Undefined value = ',undef
                0550         write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
                0551         write(6,*) ' Data valid at these times: '
                0552         ndby3 = ndatebc/3
                0553         do n = 1, ndby3*3,3
                0554          write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
                0555  1000    format(3(2x,i3,':',i8,2x,i8))
                0556         enddo
                0557         write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
3f946231fb Andr*0558         _END_MASTER( myThid )
63416ca6a5 Andr*0559        endif
e17ffcecb1 Andr*0560 
                0561        if( error ) call my_exit (101)
                0562 
63416ca6a5 Andr*0563       endif
e17ffcecb1 Andr*0564 
                0565 C---------- Read SST data if necessary -------------------------------
                0566 
598e83bd78 Andr*0567       found = .false.
                0568       nd = 2
e17ffcecb1 Andr*0569 
a6b4f97600 Jean*0570 c  If model time is not within the times of saved sst data
                0571 c  from previous call to getsst then read new data
e17ffcecb1 Andr*0572 
3f946231fb Andr*0573       timemod = dfloat(nymdmod) + dfloat(nhms)   /1000000
                0574       timebc1 = dfloat(nymdbc1) + dfloat(nhmsbc1)/1000000
                0575       timebc2 = dfloat(nymdbc2) + dfloat(nhmsbc2)/1000000
                0576 
598e83bd78 Andr*0577       if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then
                0578 
                0579        do while (.not.found .and. nd .le. ndatebc)
3f946231fb Andr*0580         timebc2 = dfloat(nymdbc(nd)) + dfloat(nhmsbc(nd))/1000000
a6b4f97600 Jean*0581         if (timebc2 .gt. timemod) then
598e83bd78 Andr*0582          nymdbc1 = nymdbc(nd-1)
                0583          nymdbc2 = nymdbc(nd)
                0584          nhmsbc1 = nhmsbc(nd-1)
                0585          nhmsbc2 = nhmsbc(nd)
819a3c957e Andr*0586          call bcdata (iunit,imbc,jmbc,nd,nd+1,sstbc1,sstbc2)
598e83bd78 Andr*0587          found = .true.
                0588         else
                0589          nd = nd + 1
                0590         endif
                0591        enddo
e17ffcecb1 Andr*0592 
a6b4f97600 Jean*0593 c  Otherwise the data from the last time in getsst surrounds the
e17ffcecb1 Andr*0594 c  current model time.
                0595 
598e83bd78 Andr*0596       else
                0597        found = .true.
                0598       endif
a6b4f97600 Jean*0599 
598e83bd78 Andr*0600       if (.not.found) then
                0601        print *, 'STOP: Could not find SST dates for model time.'
                0602        call my_finalize
                0603        call my_exit (101)
                0604       endif
e17ffcecb1 Andr*0605 
                0606 C---------- Interpolate SST data ------------------------------------
                0607 
598e83bd78 Andr*0608       call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
                0609      .                                                        fac1,fac2)
e17ffcecb1 Andr*0610 
598e83bd78 Andr*0611       do j = jm1,jm2
                0612       do i = im1,im2
a6b4f97600 Jean*0613        sst(i,j,bi,bj) = sstbc1(i+bislot,j+bjslot)*fac1
                0614      .                + sstbc2(i+bislot,j+bjslot)*fac2
598e83bd78 Andr*0615       enddo
                0616       enddo
a6b4f97600 Jean*0617 
0e97a15c40 Andr*0618 
e17ffcecb1 Andr*0619       return
                0620       end
                0621 
819a3c957e Andr*0622       subroutine bcdata (iunit,im,jm,nrec1,nrec2,field1,field2)
e17ffcecb1 Andr*0623 C************************************************************************
                0624 C
598e83bd78 Andr*0625 C!ROUTINE:      BCDATA
                0626 C!DESCRIPTION:  BCDATA reads the data from the file assigned to the
                0627 C!              passed unit number and returns data from the two times
a6b4f97600 Jean*0628 C!              surrounding the current model time.  The two record
598e83bd78 Andr*0629 C!              postitions are not assumed to be next to each other.
e17ffcecb1 Andr*0630 C
                0631 C!INPUT PARAMETERS:
598e83bd78 Andr*0632 C!      im      number of x points
                0633 C!      im      number of x points
                0634 C!      nrec1   record number of the time before the model time
                0635 C!      nrec2   record number of the time after the model time
e17ffcecb1 Andr*0636 C
                0637 C!OUTPUT PARAMETERS:
819a3c957e Andr*0638 C!      field1(im,jm)  data field before the model time
                0639 C!      field2(im,jm)  data field after the model time
e17ffcecb1 Andr*0640 C
                0641 C--------------------------------------------------------------------------
                0642       implicit none
                0643 
819a3c957e Andr*0644       integer iunit,im,jm,nrec1,nrec2
e17ffcecb1 Andr*0645 
819a3c957e Andr*0646       _RL  field1(im,jm)
                0647       _RL  field2(im,jm)
e17ffcecb1 Andr*0648 
819a3c957e Andr*0649       integer i,j
                0650       real*4 f1(im,jm), f2(im,jm)
e17ffcecb1 Andr*0651 
                0652 C--------- Read file -----------------------------------------------
                0653       read(iunit,rec=nrec1) f1
                0654       read(iunit,rec=nrec2) f2
                0655 
be3307cdae Andr*0656 #ifdef _BYTESWAPIO
819a3c957e Andr*0657       call MDS_BYTESWAPR4( im*jm, f1)
                0658       call MDS_BYTESWAPR4( im*jm, f2)
be3307cdae Andr*0659 #endif
e17ffcecb1 Andr*0660       do j=1,jm
                0661       do i=1,im
819a3c957e Andr*0662        field1(i,j) = f1(i,j)
                0663        field2(i,j) = f2(i,j)
e17ffcecb1 Andr*0664       enddo
                0665       enddo
                0666 
                0667       return
                0668       end
                0669       subroutine bcheader (iunit, ndmax, nrec,
a6b4f97600 Jean*0670      .           cname, cdscrip, im, jm, lat0, lon0, ndatebc,
598e83bd78 Andr*0671      .           nymdbc, nhmsbc, undef, error)
e17ffcecb1 Andr*0672 C************************************************************************
                0673 C
598e83bd78 Andr*0674 C!ROUTINE:     BCHEADER
a6b4f97600 Jean*0675 C!DESCRIPTION: BCHEADER reads the header from a file and returns the info.
e17ffcecb1 Andr*0676 C
                0677 C!INPUT PARAMETERS:
598e83bd78 Andr*0678 C!      iunit    unit number assigned to the data file
                0679 C!      ndmax    maximum number of date/times of the data
21eb3bf70f Jean*0680 C!      nrec     record number of the header info (or assume 1 ?)
e17ffcecb1 Andr*0681 C
                0682 C!OUTPUT PARAMETERS:
598e83bd78 Andr*0683 C!      cname         name of the data in the file header
a6b4f97600 Jean*0684 C!      cdscrip       description of the data in the file header
598e83bd78 Andr*0685 C!      im            number of x points
                0686 C!      jm            number of y points
                0687 C!      lat0          starting latitude for the data grid
                0688 C!      lon0          starting longitude for the data grid
a6b4f97600 Jean*0689 C!      ndatebc       number of date/times of the data in the file
598e83bd78 Andr*0690 C!      nymdbc(ndmax) array of dates for the data including century
                0691 C!      nhmsbc(ndmax) array of times for the data
                0692 C!      undef         value for undefined values in the data
                0693 C!      error         logical TRUE if dataset problem
e17ffcecb1 Andr*0694 C
                0695 C--------------------------------------------------------------------------
                0696       implicit none
                0697 
598e83bd78 Andr*0698       integer iunit, ndmax, nrec
                0699 
                0700       character*8  cname
                0701       character*80 cdscrip
819a3c957e Andr*0702       character*112 dummy112
                0703       integer im,jm,ndatebc,nymdbc(ndmax),nhmsbc(ndmax)
a456aa407c Andr*0704       _RL lat0,lon0,undef
598e83bd78 Andr*0705       logical error
                0706 
4f001d3c64 Andr*0707       integer i
819a3c957e Andr*0708       integer*4 im_32,jm_32
613fa3996d Andr*0709       integer*4 ndatebc_32,nhmsbc_32(ndmax),nymdbc_32(ndmax)
598e83bd78 Andr*0710       real*4 lat0_32,lon0_32,undef_32
e17ffcecb1 Andr*0711 
                0712 C--------- Read file -----------------------------------------------
                0713 
a6b4f97600 Jean*0714       read(iunit,rec=nrec,err=500) cname, cdscrip,
                0715      .     im_32, jm_32, lat0_32, lon0_32,
819a3c957e Andr*0716      .     ndatebc_32, undef_32
                0717 
be3307cdae Andr*0718 #ifdef _BYTESWAPIO
a943f0fc03 Andr*0719       call MDS_BYTESWAPI4( 1, im_32)
                0720       call MDS_BYTESWAPI4( 1, jm_32)
be3307cdae Andr*0721       call MDS_BYTESWAPR4( 1, lat0_32)
                0722       call MDS_BYTESWAPR4( 1, lon0_32)
819a3c957e Andr*0723       call MDS_BYTESWAPI4( 1, ndatebc_32)
be3307cdae Andr*0724       call MDS_BYTESWAPR4( 1, undef_32)
                0725 #endif
598e83bd78 Andr*0726 
819a3c957e Andr*0727       read(iunit,rec=nrec,err=500) dummy112,
                0728      .     (nymdbc_32(i), nhmsbc_32(i), i=1,ndatebc_32)
                0729 
598e83bd78 Andr*0730       im = im_32
                0731       jm = jm_32
                0732       lat0 = lat0_32
                0733       lon0 = lon0_32
e17ffcecb1 Andr*0734       undef = undef_32
                0735 
                0736       ndatebc = ndatebc_32
                0737       do i=1,ndatebc
819a3c957e Andr*0738 #ifdef _BYTESWAPIO
                0739       call MDS_BYTESWAPI4( 1, nymdbc_32(i))
                0740       call MDS_BYTESWAPI4( 1, nhmsbc_32(i))
                0741 #endif
e17ffcecb1 Andr*0742       nymdbc(i) = nymdbc_32(i)
                0743       nhmsbc(i) = nhmsbc_32(i)
                0744       enddo
                0745 
                0746       return
                0747   500 continue
                0748       print *, 'Error reading boundary condition from unit ',iunit
                0749       error = .true.
                0750       return
                0751       end