/[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.3 by molod, Tue Jun 8 16:42:54 2004 UTC revision 1.4 by molod, Tue Jun 8 22:26:08 2004 UTC
# Line 20  c--------------------------------------- Line 20  c---------------------------------------
20    
21         integer myTime, myIter, myThid         integer myTime, myIter, myThid
22    
23         integer i, j, L, bi, bj         integer i, j, bi, bj, biglobal, bjglobal
24         integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2         integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
25           integer nSxglobal, nSyglobal
        im1 = 1-OLx  
        im2 = sNx+OLx  
        jm1 = 1-OLy  
        jm2 = sNy+OLy  
        idim1 = 1  
        idim2 = sNx  
        jdim1 = 1  
        jdim2 = sNy  
   
26         integer ksst,kice         integer ksst,kice
   
 c Declare Local Variables  
 c -----------------------  
27         real        sstmin         real        sstmin
28         parameter ( sstmin = 273.16 )         parameter ( sstmin = 273.16 )
29    
30         integer i,j,im,jm         idim1 = 1-OLx
31           idim2 = sNx+OLx
32           jdim1 = 1-OLy
33           jdim2 = sNy+OLy
34           im1 = 1
35           im2 = sNx
36           jm1 = 1
37           jm2 = sNy
38           nSxglobal = nSx*nPx
39           nSyglobal = nSy*nPy
40    
41           call mdsfindunit( ksst, myThid )
42           call mdsfindunit( kice, myThid )
43    
44    C***********************************************************************
45    
46  C*********************************************************************         DO BJ = myByLo(myThid),myByHi(myThid)
47  C****             Interpolate Data to Current Time                ****         DO BI = myBxLo(myThid),myBxHi(myThid)
 C*********************************************************************  
48    
49         call getsst (ksst,nymd,nhms,sst,im,jm)         biglobal=bi+(myXGlobalLo-1)/im2
50         call getsice (kice,nymd,nhms,sice,im,jm)         bjglobal=bj+(myYGlobalLo-1)/jm2
51    
52           call getsst(ksst,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSx,
53         .  nSy,nSxglobal,nSyglobal,bi,bj,biglobal,bjglobal,nymd,nhms,sst)
54           call getsice(kice,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSx,
55         .  nSy,nSxglobal,nSyglobal,bi,bj,biglobal,bjglobal,nymd,nhms,sice)
56    
57  c Check for Minimum Open-Water SST  c Check for Minimum Open-Water SST
58  c --------------------------------  c --------------------------------
59         do j=1,jm         do j=jm1,jm2
60         do i=1,im         do i=im1,im2
61         if(sice(i,j).eq.0.0 .and. sst(i,j).lt.sstmin)sst(i,j) = sstmin         if(sice(i,j,bi,bj).eq.0.0 .and. sst(i,j,bi,bj).lt.sstmin)
62         .                                          sst(i,j,bi,bj) = sstmin
63         enddo         enddo
64         enddo         enddo
65    
66           ENDDO
67           ENDDO
68    
69         return         return
70         end         end
71    
72         subroutine getsice ( iunit,nymd,nhms,sice,im,jm )         subroutine getsice(iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,
73  C************************************************************************       .            nSx,nSy,nPx,nPy,bi,bj,biglobal,bjglobal,nymd,nhms,sice)
 C!GETSICE  
74  C************************************************************************  C************************************************************************
75  C  C
76  C!ROUTINE:      GETSICE  C!ROUTINE:      GETSICE
 C!PROGRAMMER:   Sharon Nebuda  
 C!DATE CODED:   May 8, 1996  
77  C!DESCRIPTION:  GETSICE returns the sea ice depth.  C!DESCRIPTION:  GETSICE returns the sea ice depth.
78  C!              This routine is adaptable for any frequency  C!              This routine is adaptable for any frequency
79  C!              data upto a daily frequency.    C!              data upto a daily frequency.  
80  C!              note: for diurnal data ndmax should be increased.  C!              note: for diurnal data ndmax should be increased.
81  C  C
82  C!INPUT PARAMETERS:  C!INPUT PARAMETERS:
83  C!              iunit   Unit number assigned to the sice data file  C!      iunit     Unit number assigned to the sice data file
84  C!              nymd    YYMMDD of the current model timestep  C!      idim1     Start dimension in x-direction
85  C!              nhms    HHMMSS of the model time  C!      idim2     End dimension in x-direction
86  C!            im     Number of x points  C!      jdim1     Start dimension in y-direction
87  C!            jm     Number of y points  C!      jdim2     End dimension in y-direction
88  C!       lattice     Grid Decomposition defined by Dynamics  C!      im1       Begin of x-direction span for filling sice
89    C!      im2       End of x-direction span for filling sice
90    C!      jm1       Begin of y-direction span for filling sice
91    C!      jm2       End of y-direction span for filling sice
92    C!      nSx       Number of processors in x-direction (local processor)
93    C!      nSy       Number of processors in y-direction (local processor)
94    C!      nPx       Number of processors in x-direction (global)
95    C!      nPy       Number of processors in y-direction (global)
96    C!      bi        Processor number in x-direction (local to processor)
97    C!      bj        Processor number in y-direction (local to processor)
98    C!      biglobal  Processor number in x-direction (global)
99    C!      bjglobal  Processor number in y-direction (global)
100    C!      nymd      YYMMDD of the current model timestep
101    C!      nhms      HHMMSS of the model time
102  C  C
103  C!OUTPUT PARAMETERS:  C!OUTPUT PARAMETERS:
104  C!              sice(im,jm)     Sea ice depth in meters  C!      sice(idim1:idim2,jdim1:jdim2,nSx,nSy) Sea ice depth in meters
105  C  C
106  C!ROUTINES CALLED:  C!ROUTINES CALLED:
107  C  C
108  C!      bcdata          Reads the data for a given unit number  C!      bcdata       Reads the data for a given unit number
109  C!      bcheader        Reads the header info for a given unit number  C!      bcheader     Reads the header info for a given unit number
110  C!     interp_time   Returns weights for linear interpolation  C!      interp_time  Returns weights for linear interpolation
111  C  C
112  C--------------------------------------------------------------------------  C--------------------------------------------------------------------------
113    
       use dynamics_lattice_module  
114        implicit none        implicit none
115        type ( dynamics_lattice_type ) lattice  #include "CPP_OPTIONS.h"
116    
117          integer iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSx,nSy
118          integer bi,bj,biglobal.bjglobal,nymd,nhms
119    
120          _RL sice(idim1:idim2,jdim1:jdim2,nSx,nSy)
121    
122  c MPI Utilities  c MPI Utilities
123  c -------------  c -------------
124        include 'mpif.h'        include 'mpif.h'
125        integer  ierror        integer  ierror
126    
127  C  Parameter statements  C Maximum number of dates in one year for the data
128        integer   ndmax   ! Maximum number of dates in one year for the data        integer   ndmax
129        parameter (ndmax = 370)        parameter (ndmax = 370)
130    
131  C  Variables passed to the routine:        character*8  cname
132        integer   iunit   ! Unit number assigned to the sea ice data file        character*80 cdscrip
133        integer   nymd    ! YYMMDD of the current model timestep        real fac1, fac2, lat0, lon0, timebc1, timebc2, timemod, undef
134        integer   nhms    ! HHMMSS of the model time        logical first, found, error
135        integer        im     ! Number of x points        integer i,j,n,nn,iyear,iyearbc,nd,ndby3,imbc,jmbc,ndatebc
136        integer        jm     ! Number of y points        integer nhmsbc(ndmax), nhmsbc1, nhmsbc2,nrec
137          integer nymdbc(ndmax),nymdbc1,nymdbc2,nymdmod
 C  Variables returned by the routine:  
       real      sice(im,jm)     ! Sea ice depth in meters  
   
 C  Variables unique to the routine:  
       character*8  cname        ! Name of the data in the file header  
       character*80 cdscrip      ! Description of the data in the file header  
       real      fac1            ! Weighted value (fraction) of the data  
                                 ! before the model time  
       real      fac2            ! Weighted value (fraction) of the data  
                                 ! after the model time  
       logical   first           ! True for first time using the dates for the  
                                 ! BC data file.  Then read in the header file.  
       logical   found           ! If false, then the data surrounding the model  
                                 ! time was not found  
       logical   error           ! TRUE if problem with data  
       integer   i,j,n,nn        ! DO loop counters  
       integer   iyear           ! Year of model  
       integer   iyearbc         ! Year of boundary condition data  
       real      lat0            ! Starting lat of the bc data set (future use)  
       real      lon0            ! Starting lon of the bc data set (future use)  
       integer   nd              ! Counter for record number of data to read  
       integer   ndby3           ! int(ndatebc/3) used for write statement  
       integer   imbc            ! IM read from the BC data  
       integer   jmbc            ! JM read from the BC data  
       integer   ndatebc         ! Number of dates in the BC file  
       integer   nhmsbc(ndmax)   ! HHMMSS of the data time (not needed currently)  
       integer   nhmsbc1 ! HHMMSS of the earlier data kept from last timestep  
       integer   nhmsbc2 ! HHMMSS of the later data kept from last timestep  
       integer   nrec            ! Record number of the header (set to 1)  
       integer   nymdbc(ndmax)   ! YYYYMMDD of each data  
       integer   nymdbc1 ! YYYYMMDD of the earlier data kept from last timestep  
       integer   nymdbc2 ! YYYYMMDD of the later data kept from last timestep  
       integer   nymdmod         ! YYYYMMDD of the current model timestep  
       real      timebc1         ! YYYYMMDD.HHMMSS of the earlier bc data  
       real      timebc2         ! YYYYMMDD.HHMMSS of the later bc data  
       real      timemod         ! YYYYMMDD.HHMMSS of the current timestep  
       real      undef           ! Undefined value for missing data  
138    
139  C Define Allocatable Arrays  C Define Allocatable Arrays
140        real, allocatable, save :: sicebc1(:,:)    ! Sea ice 0=no 1=yes from the bc data        real sicebc1(im2,jm2,nPx,nPy),sicebc2(im2,jm2,nPx,nPy)
                                                  ! of the date before the model time  
       real, allocatable, save :: sicebc2(:,:)    ! Sea ice 0=no 1=yes from bc data  
                                                  ! of the date after the model time  
141    
142  C--------- Variable Initialization ---------------------------------  C--------- Variable Initialization ---------------------------------
143    
# Line 164  c  save header info Line 148  c  save header info
148        save imbc, jmbc, lat0, lon0, ndatebc, undef, nymdbc, nhmsbc        save imbc, jmbc, lat0, lon0, ndatebc, undef, nymdbc, nhmsbc
149        save nymdbc1, nymdbc2, nhmsbc1, nhmsbc2        save nymdbc1, nymdbc2, nhmsbc1, nhmsbc2
150        save first        save first
151          save sicebc1, sicebc2
152    
153  c  this only works for between 1950-2050  c  this only works for between 1950-2050
154        if (nymd .lt. 500101) then        if (nymd .lt. 500101) then
# Line 179  c  these vaules make the iyear .ne. iyea Line 164  c  these vaules make the iyear .ne. iyea
164  c  for the first time so first isnt checked below.  c  for the first time so first isnt checked below.
165    
166        if (first) then        if (first) then
       allocate ( sicebc1(lattice%imglobal,lattice%jmglobal) )    ! Allocate Memory for sicebc1  
       allocate ( sicebc2(lattice%imglobal,lattice%jmglobal) )    ! Allocate Memory for sicebc2  
167          nymdbc(2) = 0          nymdbc(2) = 0
168          nymdbc1   = 0          nymdbc1   = 0
169          nymdbc2   = 0          nymdbc2   = 0
# Line 195  C---------- Read in Header file -------- Line 178  C---------- Read in Header file --------
178        iyearbc = nymdbc(2)/10000        iyearbc = nymdbc(2)/10000
179    
180        if( iyear.ne.iyearbc ) then        if( iyear.ne.iyearbc ) then
         if( lattice%myid.eq.0 ) then  
181    
182          close(iunit)          close(iunit)
183          open (iunit, form='unformatted', access='direct',recl=lattice%imglobal*lattice%jmglobal*4)          open (iunit, form='unformatted', access='direct',recl=im2*jm2*4)
184          nrec = 1          nrec = 1
185          call bcheader (iunit, ndmax, nrec,          call bcheader (iunit, ndmax, nrec,
186       .                 cname, cdscrip, imbc, jmbc, lat0, lon0,       .                 cname, cdscrip, imbc, jmbc, lat0, lon0,
# Line 208  C--------- Check data for Compatibility Line 190  C--------- Check data for Compatibility
190    
191  C Check for correct data in boundary condition file  C Check for correct data in boundary condition file
192         if (.not.error .and. cname.ne.'SICE') then         if (.not.error .and. cname.ne.'SICE') then
193              write(6,*) 'Wrong data in SICE boundary condition file => ',cname          write(6,*)'Wrong data in SICE boundary condition file => ',cname
194              error = .true.          error = .true.
195         endif         endif
196    
197  C Check Horizontal Resolution  C Check Horizontal Resolution
198         if (.not.error .and. imbc*jmbc.ne.lattice%imglobal*lattice%jmglobal) then         if (.not.error .and. imbc*jmbc.ne.im2*jm2) then
199              write(6,*) 'Boundary Condition Resolution DOES NOT match Model Resolution!'          write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
200              write(6,*) 'Boundary Condition Resolution:  ',imbc*jmbc          write(6,*) ' B.C. Resolution:  ',imbc*jmbc
201              write(6,*) '             Model Resolution:  ',lattice%imglobal*lattice%jmglobal          write(6,*) 'Model Resolution:  ',im2*jm2
202              error = .true.          error = .true.
203         endif         endif
204    
205  C Check Year  C Check Year
206         iyearbc = nymdbc(2)/10000         iyearbc = nymdbc(2)/10000
207         if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then         if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
208              write(6,*) 'Boundary Condition Year DOES NOT match REQUESTED Year!'          write(6,*) '     B.C. Year DOES NOT match REQUESTED Year!'
209              write(6,*) 'Boundary Condition Year:  ', iyearbc          write(6,*) '     B.C. Year:  ', iyearbc
210              write(6,*) '         Requested Year:  ', iyear          write(6,*) 'Requested Year:  ', iyear
211              error = .true.          error = .true.
212         endif         endif
213    
214         if (.not.error)   then                   if (.not.error)   then
215  C if climatology, fill dates for data with current model year  C if climatology, fill dates for data with current model year
216              if (iyearbc.eq.0) then                        if (iyearbc.eq.0) then          
217              write(6,*)           write(6,*)
218              write(6,*) 'Climatological Dataset is being used.'             write(6,*) 'Climatological Dataset is being used.'  
219              write(6,*) 'Current model year will be used to fill Header Dates'           write(6,*) 'Current model year to be used to fill Header Dates'
220                do n = 2, ndatebc-1           do n = 2, ndatebc-1
221                 nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000            nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
222                enddo           enddo
223  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
224                n = 1           n = 1
225                nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000           nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
226  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
227                n = ndatebc           n = ndatebc
228                nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000           nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
229              endif          endif
230    
231  C  Write out header info  C  Write out header info
232         write(6,*) ' Updated boundary condition data'          write(6,*) ' Updated boundary condition data'
233         write(6,*) ' ---------------------------------'          write(6,*) ' ---------------------------------'
234         write(6,*) ' Variable: ',cname          write(6,*) ' Variable: ',cname
235         write(6,*) ' Description: ',cdscrip          write(6,*) ' Description: ',cdscrip
236         write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,' Undefined value = ',undef          write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
237         write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0       .                                       ' Undefined value = ',undef
238         write(6,*) ' Data valid at these times: '          write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
239         ndby3 = ndatebc/3          write(6,*) ' Data valid at these times: '
240         do n = 1, ndby3*3,3          ndby3 = ndatebc/3
241          write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)          do n = 1, ndby3*3,3
242   1000   format(3(2x,i3,':',i8,2x,i8))           write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
243         enddo   1000    format(3(2x,i3,':',i8,2x,i8))
244         write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)          enddo
245        endif  ! End error  Check          write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
246           endif  
       endif  ! End MYID=0 Check  
   
 c Broadcast information to other PEs  
 c ----------------------------------  
 #if (mpi)  
        call mpi_bcast ( error,1,mpi_logical,0,lattice%comm,ierror )  
 #endif  
        if( error ) call my_exit (101)  
   
 #if (mpi)  
        call mpi_bcast ( ndatebc,1     ,mpi_integer,0,lattice%comm,ierror )  
        call mpi_bcast ( nymdbc,ndatebc,mpi_integer,0,lattice%comm,ierror )  
        call mpi_bcast ( nhmsbc,ndatebc,mpi_integer,0,lattice%comm,ierror )  
 #endif  
247    
248        endif  ! New Year Info Check        endif
249    
250  C---------- Read sice data if necessary -------------------------------  C---------- Read sice data if necessary -------------------------------
251    
252          found = .false.        found = .false.
253          nd = 2        nd = 2
254    
255  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
256  c  from previous call to getsice then read new data  c  from previous call to getsice then read new data
257    
258          timemod = float(nymdmod) + float(nhms)   /1000000        timemod = float(nymdmod) + float(nhms)   /1000000
259          timebc1 = float(nymdbc1) + float(nhmsbc1)/1000000        timebc1 = float(nymdbc1) + float(nhmsbc1)/1000000
260          timebc2 = float(nymdbc2) + float(nhmsbc2)/1000000        timebc2 = float(nymdbc2) + float(nhmsbc2)/1000000
261          if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then  
262          if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then
263          do while (.not.found .and. nd .le. ndatebc)  
264            timebc2 = float(nymdbc(nd)) + float(nhmsbc(nd))/1000000         do while (.not.found .and. nd .le. ndatebc)
265            if (timebc2 .gt. timemod) then                  timebc2 = float(nymdbc(nd)) + float(nhmsbc(nd))/1000000
266              nymdbc1 = nymdbc(nd-1)          if (timebc2 .gt. timemod) then  
267              nymdbc2 = nymdbc(nd)            nymdbc1 = nymdbc(nd-1)
268              nhmsbc1 = nhmsbc(nd-1)            nymdbc2 = nymdbc(nd)
269              nhmsbc2 = nhmsbc(nd)            nhmsbc1 = nhmsbc(nd-1)
270              if ( lattice%myid.eq.0 ) call bcdata (iunit, imbc, jmbc, nd, nd+1, sicebc1, sicebc2)            nhmsbc2 = nhmsbc(nd)
271  #if (mpi)            call bcdata (iunit, imbc, jmbc, nd, nd+1, sicebc1, sicebc2)
272              call mpi_bcast ( sicebc1,lattice%imglobal*lattice%jmglobal,mpi_double_precision,0,lattice%comm,ierror )            found = .true.
273              call mpi_bcast ( sicebc2,lattice%imglobal*lattice%jmglobal,mpi_double_precision,0,lattice%comm,ierror )          else
274  #endif            nd = nd + 1
275              found = .true.          endif
276            else         enddo
             nd = nd + 1  
           endif  
         enddo  
277    
278  c  Otherwise the data from the last time in getsice surrounds the  c  Otherwise the data from the last time in getsice surrounds the
279  c  current model time.  c  current model time.
280    
281          else        else
282            found = .true.         found = .true.
283          endif        endif
284    
285            if (.not.found) then        if (.not.found) then
286              if( lattice%myid.eq.0 ) print *, 'STOP: Could not find SICE boundary condition dates surrounding the model time.'         print *, 'STOP: Could not find SICE dates for model time.'
287              call my_finalize         call my_finalize
288              call my_exit (101)         call my_exit (101)
289            endif        endif
290    
291  C---------- Interpolate sice data ------------------------------------  C---------- Interpolate sice data ------------------------------------
292    
293            call interp_time ( nymdmod,nhms, nymdbc1,nhmsbc1, nymdbc2,nhmsbc2, fac1,fac2 )        call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
294         .                                                       fac1,fac2)
295    
296            do j = 1, jm        do j = jm1,jm2
297            do i = 1, im        do i = im1,im2
298              sice(i,j) = sicebc1( lattice%iglobal(i),lattice%jglobal(j) )*fac1         sice(i,j,bi,bj) = sicebc1(i,j,biglobal,bjglobal)*fac1
299       .                + sicebc2( lattice%iglobal(i),lattice%jglobal(j) )*fac2       .                 + sicebc2(i,j,biglobal,bjglobal)*fac2
300  c average to 0 or 1  c average to 0 or 1
301  c -----------------  c -----------------
302              if (sice(i,j) .ge. 0.5) then         if (sice(i,j,bi,bj) .ge. 0.5) then
303                  sice(i,j) = 1.          sice(i,j,bi,bj) = 1.
304              else         else
305                  sice(i,j) = 0.          sice(i,j,bi,bj) = 0.
306              endif         endif
307            enddo        enddo
308            enddo        enddo
309    
310  C---------- Fill sice with depth of ice ------------------------------------  C---------- Fill sice with depth of ice ------------------------------------
311          do j = jm1,jm2
312            do j = 1, jm        do i = im1,im2
313            do i = 1, im         if (sice(i,j,bi,bj) .eq. 1.) then                
314              if (sice(i,j) .eq. 1.) then         ! sea ice present          sice(i,j,bi,bj) = 3.
315                  sice(i,j) = 3.         endif
316              endif        enddo
317            enddo        enddo
           enddo  
   
318  C---------------------------------------------------------------------------  C---------------------------------------------------------------------------
319    
320        return        return
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22