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

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

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

revision 1.3 by molod, Mon Jun 7 18:04:06 2004 UTC revision 1.21 by molod, Mon Mar 7 20:01:01 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4        subroutine fizhi_init_veg(mythid,vegdata,im,jm,maxbi,maxbj,maxtyp,  #include "FIZHI_OPTIONS.h"
5       . surftype,tilefrac,igrd,ityp,chfr )  
6          subroutine fizhi_init_veg(mythid,vegdata,im,jm,Nsx,Nsy,Nxg,Nyg,
7         . maxtyp,nchp,nchptot,nchpland,lons,lats,surftype,tilefrac,
8         . igrd,ityp,chfr,chlt,chlon)
9  C***********************************************************************  C***********************************************************************
10  C Subroutine fizhi_init_veg - routine to read in the land surface types,  C Subroutine fizhi_init_veg - routine to read in the land surface types,
11  C      interpolate to the models grid, and set up tile space for use by  C      interpolate to the models grid, and set up tile space for use by
# Line 13  C INPUT: Line 16  C INPUT:
16  C  C
17  C mythid   - thread number (processor number)  C mythid   - thread number (processor number)
18  C vegdata  - Character*40 Vegetation Dataset name  C vegdata  - Character*40 Vegetation Dataset name
19  C im       - model grid longitude dimension  C im       - longitude dimension
20  C jm       - model grid latitude dimension (number of lat. points)  C jm       - latitude dimension (number of lat. points)
21  C maxbi    - Number of processors in x-direction  C Nsx      - Number of processors in x-direction
22  C maxbj    - Number of processors in y-direction  C Nsy      - Number of processors in y-direction
23  C maxtyp   - maximum allowable number of land surface types per grid box  C maxtyp   - maximum allowable number of land surface types per grid box
24  C nchpmax  - integer maximum per-processor number of tiles in tile space  C nchp     - integer per-processor number of tiles in tile space
25    C lons     - longitude in degrees [im,jm,nSx,nSy]
26    C lats     - latitude in degrees [im,jm,nSx,nSy]
27  C  C
28  C OUTPUT:  C OUTPUT:
29  C  C
30  C surftype - integer array of land surface types [im,jm,maxtyp,bi,bj]  C surftype - integer array of land surface types [im,jm,maxtyp,Nsx,Nsy]
31  C tilefrac - real array of corresponding land surface type fractions  C tilefrac - real array of corresponding land surface type fractions
32  C            [im,jm,maxtyp,bi,bj]  C            [im,jm,maxtyp,Nsx,Nsy]
33  C igrd     - integer array in tile space of grid point number for each  C igrd     - integer array in tile space of grid point number for each
34  C            tile [nchp,bi,bj]  C            tile [nchp,Nsx,Nsy]
35  C ityp     - integer array in tile space of land surface type for each  C ityp     - integer array in tile space of land surface type for each
36  C            tile [nchp,bi,bj]  C            tile [nchp,Nsx,Nsy]
37  C chfr     - real array in tile space of land surface type fraction for  C chfr     - real array in tile space of land surface type fraction for
38  C            each tile [nchp,bi,bj]  C            each tile [nchp,Nsx,Nsy]
39  C  C
40  C NOTES:  C NOTES:
41  C       Vegetation type as follows:  C       Vegetation type as follows:
# Line 49  C*************************************** Line 54  C***************************************
54        implicit none        implicit none
55  #include "EEPARAMS.h"  #include "EEPARAMS.h"
56    
57        integer mythid,im,jm,maxtyp,nchpmax,maxbi,maxbj        integer mythid,im,jm,maxtyp,nchp,Nsx,Nsy,Nxg,Nyg
58        integer surftype(im,jm,maxtyp,maxbi,maxbj)        integer nchptot(Nsx,Nsy), nchpland(Nsx,Nsy)
59        integer igrd(nchpmax,bi,bj),ityp(nchpmax,maxbi,maxbj)        integer surftype(im,jm,maxtyp,Nsx,Nsy)
60        real tilefrac(im,jm,maxtyp,maxbi,maxbj)        integer igrd(nchp,Nsx,Nsy),ityp(nchp,Nsx,Nsy)
61        real chfr(nchpmax,maxbi,maxbj)        _RL tilefrac(im,jm,maxtyp,Nsx,Nsy)
62          _RL lats(im,jm,nSx,nSy), lons(im,jm,nSx,nSy)
63          _RL chfr(nchp,Nsx,Nsy),chlt(nchp,Nsx,Nsy),chlon(nchp,Nsx,Nsy)
64        character*40 vegdata        character*40 vegdata
65        integer imdata,jmdata,bidata,bjdata        integer imdata,jmdata,Nxgdata,Nygdata
66        integer nchp,nchpland        integer biglobal,bjglobal
67    
68        integer*4 im_32, jm_32, bi_32, bj_32        integer*4 im_32, jm_32, Nxg_32, Nyg_32
69        integer*4 iveg_32(im,jm,maxtyp,bi,bj)        integer*4 iveg_32(im,jm,maxtyp,Nxg,Nyg)
70           real*4  veg_32(im,jm,maxtyp,bi,bj)           real*4  veg_32(im,jm,maxtyp,Nxg,Nyg)
71    
72        integer i,j,k,bi,bj,ierr1,kveg        integer i,j,k,bi,bj,ierr1,kveg
73    
# Line 69  C*************************************** Line 76  C***************************************
76        open(kveg,file=vegdata,form='unformatted',access='sequential',        open(kveg,file=vegdata,form='unformatted',access='sequential',
77       .                      iostat=ierr1)       .                      iostat=ierr1)
78        if( ierr1.eq.0 ) then        if( ierr1.eq.0 ) then
79            read(kveg)im_32,jm_32,bi_32,bj_32,IVEG_32,VEG_32            read(kveg)im_32,jm_32,Nxg_32,Nyg_32,IVEG_32,VEG_32
80        else        else
81         print *         print *
82         print *, 'Veg Dataset: ',vegdata,' not found!'         print *, 'Veg Dataset: ',vegdata,' not found!'
# Line 77  C*************************************** Line 84  C***************************************
84         call exit(101)         call exit(101)
85        endif        endif
86        close(kveg)        close(kveg)
87    
88        IF (myThid.eq.1) THEN        IF (myThid.eq.1) THEN
89        imdata = im_32        imdata = im_32
90        jmdata = jm_32        jmdata = jm_32
91        bidata = bi_32        Nxgdata = Nxg_32
92        bjdata = bj_32        Nygdata = Nyg_32
93        if( (imdata.ne.im) .or. (jmdata.ne.jm) .or.        if( (imdata.ne.im) .or. (jmdata.ne.jm) .or.
94       .                        (bi.ne.bidata) .or. (bjdata.ne.bj) ) then       .                     (Nxgdata.ne.Nxg) .or. (Nygdata.ne.Nyg) ) then
95         print *         print *
96         print *, 'Veg Data Resolution is Incorrect! '         print *, 'Veg Data Resolution is Incorrect! '
97         print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata         print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata
98         print *,' Model Bij: ',bi,'x',bj,' Data Bij: ',bidata,'x',bjdata         print *,' Model Nxg Nyg: ',Nxg,' ',Nyg,' Data Nxg Nyg: ',Nxgdata,
99         .                    ' ',Nygdata
100         print *         print *
101         call exit(102)         call exit(102)
102        ENDIF        ENDIF
103                                                                                          ENDIF
       imdata = im_32  
       jmdata = jm_32  
       bidata = bi_32  
       bjdata = bj_32  
104    
105        DO BJ = myByLo(myThid), myByHi(myThid)        DO BJ = myByLo(myThid), myByHi(myThid)
106        DO BI = myBxLo(myThid), myBxHi(myThid)        DO BI = myBxLo(myThid), myBxHi(myThid)
107                                                                                    
108          biglobal=bi+(myXGlobalLo-1)/im
109          bjglobal=bj+(myYGlobalLo-1)/jm
110    #if defined( _BYTESWAPIO ) && defined( ALLOW_MDSIO )
111          call MDS_BYTESWAPR4(im*jm*maxtyp,veg_32(1,1,1,biglobal,bjglobal))
112    #endif
113        do k = 1,maxtyp        do k = 1,maxtyp
114        do j = 1,jm_32        do j = 1,jm
115        do i = 1,im_32        do i = 1,im
116         surftype(i,j,k,bi,bj) = iveg_32(i,j,k,bi,bj)         surftype(i,j,k,bi,bj) = iveg_32(i,j,k,biglobal,bjglobal)
117         tilefrac(i,j,k,bi,bj) = veg_32(i,j,k,bi,bj)         tilefrac(i,j,k,bi,bj) = veg_32(i,j,k,biglobal,bjglobal)
118        enddo        enddo
119        enddo        enddo
120        enddo        enddo
121    
122  c   create chip arrays for :        ENDDO
123          ENDDO
124    
125    c     create chip arrays for :
126  c      igrd :  grid index  c      igrd :  grid index
127  c      ityp :  veg. type  c      ityp :  veg. type
128  c      chfr :  vegetation fraction  c      chfr :  vegetation fraction
129    c      chlon:  chip longitude
130    c      chlt :  chip latitude
131    
132  c  nchplnd<=nchp is the actual number of land chips  c     nchpland<=nchptot is the actual number of land chips
   
 c land points  
 c -----------  
       nchplnd = 0  
       do k=1,maxtyp  
       do j=1,jm  
       do i=1,im  
       if(surftype(i,j,k,bi,bj).lt.100.and.  
      .                                  tilefrac(i,j,k,bi,bj).gt.0.)then  
              nchplnd  = nchplnd + 1  
        igrd (nchplnd,bi,bj) = i + (j-1)*im  
        ityp (nchplnd,bi,bj) = surftype(i,j,k,bi,bj)  
        chfr (nchplnd,bi,bj) = tilefrac(i,j,k,bi,bj)  
       endif  
       enddo  
       enddo  
       enddo  
133    
134  c ocean points        DO BJ = myByLo(myThid), myByHi(myThid)
135  c ------------          DO BI = myBxLo(myThid), myBxHi(myThid)
       nchp = nchplnd  
   
       do k=1,maxtyp  
       do j=1,jm  
       do i=1,im  
       if(surftype(i,j,k,bi,bj).ge.100 .and.  
      .                                  tilefrac(i,j,k,bi,bj).gt.0.)then  
              nchp  = nchp + 1  
        igrd (nchp,bi,bj) = i + (j-1)*im  
        ityp (nchp,bi,bj) = surftype(i,j,k,bi,bj)  
        chfr (nchp,bi,bj) = tilefrac(i,j,k,bi,bj)  
       endif  
       enddo  
       enddo  
       enddo  
136    
137        print *, 'bi ',bi,' bj ',bj  c         land points
138        print *, 'Number of Total Tiles: ',nchp  c         -----------
139        print *, 'Number of Land  Tiles: ',nchplnd            nchpland(bi,bj) = 0
140        print *            do k=1,maxtyp
141                do j=1,jm
142                  do i=1,im
143                    if(surftype(i,j,k,bi,bj).lt.100 .and.
144         .               tilefrac(i,j,k,bi,bj).gt.0.) then
145                      nchpland(bi,bj)  = nchpland(bi,bj) + 1
146                      igrd (nchpland(bi,bj),bi,bj) = i + (j-1)*im
147                      ityp (nchpland(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
148                      chfr (nchpland(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
149                      chlon(nchpland(bi,bj),bi,bj) = lons(i,j,bi,bj)
150                      chlt (nchpland(bi,bj),bi,bj) = lats(i,j,bi,bj)
151                    endif
152                  enddo
153                enddo
154              enddo
155              
156    c         ocean points
157    c         ------------
158              nchptot(bi,bj) = nchpland(bi,bj)
159              
160              do k=1,maxtyp
161                do j=1,jm
162                  do i=1,im
163                    if(surftype(i,j,k,bi,bj).ge.100 .and.
164         .               tilefrac(i,j,k,bi,bj).gt.0.) then
165                      nchptot(bi,bj)  = nchptot(bi,bj) + 1
166                      igrd (nchptot(bi,bj),bi,bj) = i + (j-1)*im
167                      ityp (nchptot(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
168                      chfr (nchptot(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
169                      chlon(nchptot(bi,bj),bi,bj) = lons(i,j,bi,bj)
170                      chlt (nchptot(bi,bj),bi,bj) = lats(i,j,bi,bj)
171                    endif
172                  enddo
173                enddo
174              enddo
175              
176              print *,'No of Total Tiles for bi=',bi,': ',nchptot(bi,bj)
177              print *,'No of Land  Tiles for bi=',bi,': ',nchpland(bi,bj)
178    
179        ENDDO          ENDDO
180        ENDDO        ENDDO
181    
182        RETURN        RETURN
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22