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

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

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


Revision 1.14 - (hide annotations) (download)
Thu Jul 22 15:56:01 2004 UTC (19 years, 11 months ago) by edhill
Branch: MAIN
Changes since 1.13: +90 -51 lines
 o first cut at using MNC to read FIZHI vegetation type data

1 edhill 1.14 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_init_veg.F,v 1.13 2004/07/20 16:24:49 molod Exp $
2 molod 1.1 C $Name: $
3    
4 edhill 1.14 #include "FIZHI_OPTIONS.h"
5    
6 molod 1.6 subroutine fizhi_init_veg(mythid,vegdata,im,jm,Nsx,Nsy,Nxg,Nyg,
7 molod 1.12 . maxtyp,nchp,nchpland,lons,lats,surftype,tilefrac,igrd,ityp,chfr,
8     . chlt,chlon)
9 molod 1.1 C***********************************************************************
10 molod 1.2 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
12     C the land surface model, the albedo calculation and the surface
13     C roughness calculation.
14 molod 1.1 C
15     C INPUT:
16     C
17     C mythid - thread number (processor number)
18 molod 1.2 C vegdata - Character*40 Vegetation Dataset name
19 molod 1.7 C im - longitude dimension
20     C jm - latitude dimension (number of lat. points)
21 molod 1.4 C Nsx - Number of processors in x-direction
22     C Nsy - Number of processors in y-direction
23 molod 1.1 C maxtyp - maximum allowable number of land surface types per grid box
24 molod 1.6 C nchp - integer per-processor number of tiles in tile space
25 molod 1.7 C lons - longitude in degrees [im,jm,nSx,nSy]
26     C lats - latitude in degrees [im,jm,nSx,nSy]
27 molod 1.1 C
28     C OUTPUT:
29     C
30 molod 1.4 C surftype - integer array of land surface types [im,jm,maxtyp,Nsx,Nsy]
31 molod 1.1 C tilefrac - real array of corresponding land surface type fractions
32 molod 1.4 C [im,jm,maxtyp,Nsx,Nsy]
33 molod 1.1 C igrd - integer array in tile space of grid point number for each
34 molod 1.4 C tile [nchp,Nsx,Nsy]
35 molod 1.1 C ityp - integer array in tile space of land surface type for each
36 molod 1.4 C tile [nchp,Nsx,Nsy]
37 molod 1.1 C chfr - real array in tile space of land surface type fraction for
38 molod 1.4 C each tile [nchp,Nsx,Nsy]
39 molod 1.1 C
40     C NOTES:
41     C Vegetation type as follows:
42     C 1: BROADLEAF EVERGREEN TREES
43     C 2: BROADLEAF DECIDUOUS TREES
44     C 3: NEEDLELEAF TREES
45     C 4: GROUND COVER
46     C 5: BROADLEAF SHRUBS
47     C 6: DWARF TREES (TUNDRA)
48     C 7: BARE SOIL
49     C 8: DESERT
50     C 9: GLACIER
51     C 10: DARK DESERT
52     C 100: OCEAN
53     C***********************************************************************
54     implicit none
55 molod 1.5 #include "EEPARAMS.h"
56 molod 1.1
57 molod 1.12 integer mythid,im,jm,maxtyp,nchp,nchpland,Nsx,Nsy,Nxg,Nyg
58 molod 1.4 integer surftype(im,jm,maxtyp,Nsx,Nsy)
59 molod 1.6 integer igrd(nchp,Nsx,Nsy),ityp(nchp,Nsx,Nsy)
60 molod 1.10 _RL tilefrac(im,jm,maxtyp,Nsx,Nsy)
61     _RL lats(im,jm,nSx,nSy), lons(im,jm,nSx,nSy)
62     _RL chfr(nchp,Nsx,Nsy),chlt(nchp,Nsx,Nsy),chlon(nchp,Nsx,Nsy)
63 molod 1.1 character*40 vegdata
64 molod 1.6 integer imdata,jmdata,Nxgdata,Nygdata
65 molod 1.12 integer nchplocal,biglobal,bjglobal
66 molod 1.1
67 molod 1.6 integer*4 im_32, jm_32, Nxg_32, Nyg_32
68     integer*4 iveg_32(im,jm,maxtyp,Nxg,Nyg)
69     real*4 veg_32(im,jm,maxtyp,Nxg,Nyg)
70 molod 1.1
71 molod 1.3 integer i,j,k,bi,bj,ierr1,kveg
72 molod 1.1
73 edhill 1.14 #ifdef ALLOW_MNC
74     character*(MAX_LEN_FNAM) fizhi_veg_bn
75     #endif
76    
77     C Allow for MDSIO format if someday needed
78     #ifdef ALLOW_MDSIO
79     IF ( .FALSE. ) THEN
80    
81 molod 1.1 call mdsfindunit( kveg, myThid )
82 molod 1.2 close(kveg)
83     open(kveg,file=vegdata,form='unformatted',access='sequential',
84     . iostat=ierr1)
85 molod 1.1 if( ierr1.eq.0 ) then
86 molod 1.6 read(kveg)im_32,jm_32,Nxg_32,Nyg_32,IVEG_32,VEG_32
87 molod 1.1 else
88     print *
89     print *, 'Veg Dataset: ',vegdata,' not found!'
90     print *
91 molod 1.2 call exit(101)
92 molod 1.1 endif
93 molod 1.2 close(kveg)
94 edhill 1.14
95 molod 1.3 IF (myThid.eq.1) THEN
96 molod 1.2 imdata = im_32
97     jmdata = jm_32
98 molod 1.6 Nxgdata = Nxg_32
99     Nygdata = Nyg_32
100 molod 1.2 if( (imdata.ne.im) .or. (jmdata.ne.jm) .or.
101 molod 1.6 . (Nxgdata.ne.Nxg) .or. (Nygdata.ne.Nyg) ) then
102 molod 1.2 print *
103     print *, 'Veg Data Resolution is Incorrect! '
104     print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata
105 molod 1.6 print *,' Model Nxg Nyg: ',Nxg,' ',Nyg,' Data Nxg Nyg: ',Nxgdata,
106     . ' ',Nygdata
107 molod 1.2 print *
108     call exit(102)
109 molod 1.1 ENDIF
110 molod 1.12 ENDIF
111    
112 molod 1.5 DO BJ = myByLo(myThid), myByHi(myThid)
113     DO BI = myBxLo(myThid), myBxHi(myThid)
114    
115     biglobal=bi+(myXGlobalLo-1)/im
116     bjglobal=bj+(myYGlobalLo-1)/jm
117 edhill 1.14 #if defined( _BYTESWAPIO ) && defined( ALLOW_MDSIO )
118 molod 1.13 call MDS_BYTESWAPR4(im*jm*maxtyp,veg_32(1,1,1,biglobal,bjglobal))
119     #endif
120 molod 1.2 do k = 1,maxtyp
121 molod 1.4 do j = 1,jm
122     do i = 1,im
123 molod 1.5 surftype(i,j,k,bi,bj) = iveg_32(i,j,k,biglobal,bjglobal)
124     tilefrac(i,j,k,bi,bj) = veg_32(i,j,k,biglobal,bjglobal)
125 molod 1.2 enddo
126 molod 1.1 enddo
127     enddo
128    
129 edhill 1.14 ENDDO
130     ENDDO
131    
132     ENDIF
133     #endif
134    
135     #ifdef ALLOW_MNC
136     _BEGIN_MASTER( myThid )
137    
138     do i = 1,MAX_LEN_FNAM
139     fizhi_veg_bn(i:i) = ' '
140     enddo
141    
142     C The following base name should be handled by some sort of input
143     C name parameter in FIZHI_READPARMS() plus a possible size.
144    
145     C Set the base name 1234567890
146     fizhi_veg_bn(1:10) = 'fizhi_veg '
147    
148     CALL MNC_CW_I_R('I', fizhi_veg_bn, 0,0,
149     & 'surftype',surftype,myThid)
150     CALL MNC_CW_RL_R('R', fizhi_veg_bn, 0,0,
151     & 'tilefrac', tilefrac, myThid)
152    
153     _END_MASTER( myThid )
154     #endif
155    
156     c create chip arrays for :
157 molod 1.1 c igrd : grid index
158     c ityp : veg. type
159     c chfr : vegetation fraction
160 molod 1.7 c chlon: chip longitude
161     c chlt : chip latitude
162 molod 1.1
163 edhill 1.14 c nchpland<=nchplocal is the actual number of land chips
164 molod 1.1
165 edhill 1.14 DO BJ = myByLo(myThid), myByHi(myThid)
166     DO BI = myBxLo(myThid), myBxHi(myThid)
167 molod 1.1
168 edhill 1.14 c land points
169     c -----------
170     nchpland = 0
171     do k=1,maxtyp
172     do j=1,jm
173     do i=1,im
174     if(surftype(i,j,k,bi,bj).lt.100 .and.
175     . tilefrac(i,j,k,bi,bj).gt.0.) then
176     nchpland = nchpland + 1
177     igrd (nchpland,bi,bj) = i + (j-1)*im
178     ityp (nchpland,bi,bj) = surftype(i,j,k,bi,bj)
179     chfr (nchpland,bi,bj) = tilefrac(i,j,k,bi,bj)
180     chlon(nchpland,bi,bj) = lons(i,j,bi,bj)
181     chlt (nchpland,bi,bj) = lats(i,j,bi,bj)
182     endif
183     enddo
184     enddo
185     enddo
186    
187     c ocean points
188     c ------------
189     nchplocal = nchpland
190    
191     do k=1,maxtyp
192     do j=1,jm
193     do i=1,im
194     if(surftype(i,j,k,bi,bj).ge.100 .and.
195     . tilefrac(i,j,k,bi,bj).gt.0.) then
196     nchplocal = nchplocal + 1
197     igrd (nchplocal,bi,bj) = i + (j-1)*im
198     ityp (nchplocal,bi,bj) = surftype(i,j,k,bi,bj)
199     chfr (nchplocal,bi,bj) = tilefrac(i,j,k,bi,bj)
200     chlon(nchplocal,bi,bj) = lons(i,j,bi,bj)
201     chlt (nchplocal,bi,bj) = lats(i,j,bi,bj)
202     endif
203     enddo
204     enddo
205     enddo
206    
207     print *, 'bi ',bi,' bj ',bj
208     print *, 'Number of Total Tiles: ',nchplocal
209     print *, 'Number of Land Tiles: ',nchpland
210     print *
211    
212     ENDDO
213 molod 1.2 ENDDO
214 edhill 1.14
215 molod 1.1 RETURN
216     END

  ViewVC Help
Powered by ViewVC 1.1.22