/[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.18 by molod, Fri Oct 22 19:20:55 2004 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 48  C                100:  OCEAN Line 53  C                100:  OCEAN
53  C***********************************************************************  C***********************************************************************
54        implicit none        implicit none
55  #include "EEPARAMS.h"  #include "EEPARAMS.h"
56    #include "PARAMS.h"
57    
58        integer mythid,im,jm,maxtyp,nchpmax,maxbi,maxbj        integer mythid,im,jm,maxtyp,nchp,Nsx,Nsy,Nxg,Nyg
59        integer surftype(im,jm,maxtyp,maxbi,maxbj)        integer nchptot(Nsx,Nsy), nchpland(Nsx,Nsy)
60        integer igrd(nchpmax,bi,bj),ityp(nchpmax,maxbi,maxbj)        integer surftype(im,jm,maxtyp,Nsx,Nsy)
61        real tilefrac(im,jm,maxtyp,maxbi,maxbj)        integer igrd(nchp,Nsx,Nsy),ityp(nchp,Nsx,Nsy)
62        real chfr(nchpmax,maxbi,maxbj)        _RL tilefrac(im,jm,maxtyp,Nsx,Nsy)
63          _RL lats(im,jm,nSx,nSy), lons(im,jm,nSx,nSy)
64          _RL chfr(nchp,Nsx,Nsy),chlt(nchp,Nsx,Nsy),chlon(nchp,Nsx,Nsy)
65        character*40 vegdata        character*40 vegdata
66        integer imdata,jmdata,bidata,bjdata        integer imdata,jmdata,Nxgdata,Nygdata
67        integer nchp,nchpland        integer biglobal,bjglobal
68    
69        integer*4 im_32, jm_32, bi_32, bj_32        integer*4 im_32, jm_32, Nxg_32, Nyg_32
70        integer*4 iveg_32(im,jm,maxtyp,bi,bj)        integer*4 iveg_32(im,jm,maxtyp,Nxg,Nyg)
71           real*4  veg_32(im,jm,maxtyp,bi,bj)           real*4  veg_32(im,jm,maxtyp,Nxg,Nyg)
72    
73        integer i,j,k,bi,bj,ierr1,kveg        integer i,j,k,bi,bj,ierr1,kveg
74    
75    #ifdef ALLOW_MNC
76          character*(MAX_LEN_FNAM) fizhi_veg_bn
77    #endif
78    
79    C     Allow for MDSIO format if someday needed
80    #ifdef ALLOW_MDSIO
81          IF ( .FALSE. ) THEN
82    
83        call mdsfindunit( kveg, myThid )        call mdsfindunit( kveg, myThid )
84        close(kveg)        close(kveg)
85        open(kveg,file=vegdata,form='unformatted',access='sequential',        open(kveg,file=vegdata,form='unformatted',access='sequential',
86       .                      iostat=ierr1)       .                      iostat=ierr1)
87        if( ierr1.eq.0 ) then        if( ierr1.eq.0 ) then
88            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
89        else        else
90         print *         print *
91         print *, 'Veg Dataset: ',vegdata,' not found!'         print *, 'Veg Dataset: ',vegdata,' not found!'
# Line 77  C*************************************** Line 93  C***************************************
93         call exit(101)         call exit(101)
94        endif        endif
95        close(kveg)        close(kveg)
96    
97        IF (myThid.eq.1) THEN        IF (myThid.eq.1) THEN
98        imdata = im_32        imdata = im_32
99        jmdata = jm_32        jmdata = jm_32
100        bidata = bi_32        Nxgdata = Nxg_32
101        bjdata = bj_32        Nygdata = Nyg_32
102        if( (imdata.ne.im) .or. (jmdata.ne.jm) .or.        if( (imdata.ne.im) .or. (jmdata.ne.jm) .or.
103       .                        (bi.ne.bidata) .or. (bjdata.ne.bj) ) then       .                     (Nxgdata.ne.Nxg) .or. (Nygdata.ne.Nyg) ) then
104         print *         print *
105         print *, 'Veg Data Resolution is Incorrect! '         print *, 'Veg Data Resolution is Incorrect! '
106         print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata         print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata
107         print *,' Model Bij: ',bi,'x',bj,' Data Bij: ',bidata,'x',bjdata         print *,' Model Nxg Nyg: ',Nxg,' ',Nyg,' Data Nxg Nyg: ',Nxgdata,
108         .                    ' ',Nygdata
109         print *         print *
110         call exit(102)         call exit(102)
111        ENDIF        ENDIF
112                                                                                          ENDIF
       imdata = im_32  
       jmdata = jm_32  
       bidata = bi_32  
       bjdata = bj_32  
113    
114        DO BJ = myByLo(myThid), myByHi(myThid)        DO BJ = myByLo(myThid), myByHi(myThid)
115        DO BI = myBxLo(myThid), myBxHi(myThid)        DO BI = myBxLo(myThid), myBxHi(myThid)
116                                                                                    
117          biglobal=bi+(myXGlobalLo-1)/im
118          bjglobal=bj+(myYGlobalLo-1)/jm
119    #if defined( _BYTESWAPIO ) && defined( ALLOW_MDSIO )
120          call MDS_BYTESWAPR4(im*jm*maxtyp,veg_32(1,1,1,biglobal,bjglobal))
121    #endif
122        do k = 1,maxtyp        do k = 1,maxtyp
123        do j = 1,jm_32        do j = 1,jm
124        do i = 1,im_32        do i = 1,im
125         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)
126         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)
127        enddo        enddo
128        enddo        enddo
129        enddo        enddo
130    
131  c   create chip arrays for :        ENDDO
132          ENDDO
133    
134          ENDIF
135    #endif
136    
137    #ifdef ALLOW_MNC
138          IF (useMNC) THEN
139            _BEGIN_MASTER( myThid )
140          
141            do i = 1,MAX_LEN_FNAM
142              fizhi_veg_bn(i:i) = ' '
143            enddo
144          
145    C       The following base name should be handled by some sort of input
146    C       name parameter in FIZHI_READPARMS() plus a possible size.
147          
148    C       Set the base name     1234567890
149            fizhi_veg_bn(1:10) = 'fizhi_veg '
150          
151            CALL MNC_CW_I_R('I', fizhi_veg_bn, 0,0,
152         &       'surftype',surftype,myThid)
153            CALL MNC_CW_RL_R('R', fizhi_veg_bn, 0,0,
154         &       'tilefrac', tilefrac, myThid)
155          
156            _END_MASTER( myThid )
157          ENDIF
158    #endif
159    
160    c     create chip arrays for :
161  c      igrd :  grid index  c      igrd :  grid index
162  c      ityp :  veg. type  c      ityp :  veg. type
163  c      chfr :  vegetation fraction  c      chfr :  vegetation fraction
164    c      chlon:  chip longitude
165    c      chlt :  chip latitude
166    
167  c  nchplnd<=nchp is the actual number of land chips  c     nchpland<=nchptot is the actual number of land chips
168    
169  c land points        DO BJ = myByLo(myThid), myByHi(myThid)
170  c -----------          DO BI = myBxLo(myThid), myBxHi(myThid)
       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  
   
 c ocean points  
 c ------------  
       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  
171    
172        print *, 'bi ',bi,' bj ',bj  c         land points
173        print *, 'Number of Total Tiles: ',nchp  c         -----------
174        print *, 'Number of Land  Tiles: ',nchplnd            nchpland(bi,bj) = 0
175        print *            do k=1,maxtyp
176                do j=1,jm
177                  do i=1,im
178                    if(surftype(i,j,k,bi,bj).lt.100 .and.
179         .               tilefrac(i,j,k,bi,bj).gt.0.) then
180                      nchpland(bi,bj)  = nchpland(bi,bj) + 1
181                      igrd (nchpland(bi,bj),bi,bj) = i + (j-1)*im
182                      ityp (nchpland(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
183                      chfr (nchpland(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
184                      chlon(nchpland(bi,bj),bi,bj) = lons(i,j,bi,bj)
185                      chlt (nchpland(bi,bj),bi,bj) = lats(i,j,bi,bj)
186                    endif
187                  enddo
188                enddo
189              enddo
190              
191    c         ocean points
192    c         ------------
193              nchptot(bi,bj) = nchpland(bi,bj)
194              
195              do k=1,maxtyp
196                do j=1,jm
197                  do i=1,im
198                    if(surftype(i,j,k,bi,bj).ge.100 .and.
199         .               tilefrac(i,j,k,bi,bj).gt.0.) then
200                      nchptot(bi,bj)  = nchptot(bi,bj) + 1
201                      igrd (nchptot(bi,bj),bi,bj) = i + (j-1)*im
202                      ityp (nchptot(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
203                      chfr (nchptot(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
204                      chlon(nchptot(bi,bj),bi,bj) = lons(i,j,bi,bj)
205                      chlt (nchptot(bi,bj),bi,bj) = lats(i,j,bi,bj)
206                    endif
207                  enddo
208                enddo
209              enddo
210              
211              print *,'Number of Total Tiles for bi=',bi,': ',nchptot(bi,bj)
212              print *,'Number of Land  Tiles for bi=',bi,': ',nchpland(bi,bj)
213    
214            ENDDO
215        ENDDO        ENDDO
216        ENDDO        
   
217        RETURN        RETURN
218        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22