pkg/fizhi/fizhi_init_veg.F

f9c065a827 Ed H*0001 #include "FIZHI_OPTIONS.h"
                0002
d58f306b82 Jean*0003       SUBROUTINE FIZHI_INIT_VEG(myThid,vegdata,im,jm,nSx,nSy,Nxg,Nyg,
                0004      & maxtyp,nchp,nchptot,nchpland,lons,lats,surftype,tilefrac,
                0005      & igrd,ityp,chfr,chlt,chlon)
614dbd872e Andr*0006 C***********************************************************************
45815db68b Andr*0007 C Subroutine fizhi_init_veg - routine to read in the land surface types,
d58f306b82 Jean*0008 C      interpolate to the models grid, and set up tile space for use by
                0009 C      the land surface model, the albedo calculation and the surface
45815db68b Andr*0010 C      roughness calculation.
614dbd872e Andr*0011 C
                0012 C INPUT:
d58f306b82 Jean*0013 C
                0014 C myThid   - thread number (processor number)
45815db68b Andr*0015 C vegdata  - Character*40 Vegetation Dataset name
dbae14396f Andr*0016 C im       - longitude dimension
                0017 C jm       - latitude dimension (number of lat. points)
d58f306b82 Jean*0018 C nSx      - Number of processors in x-direction
                0019 C nSy      - Number of processors in y-direction
614dbd872e Andr*0020 C maxtyp   - maximum allowable number of land surface types per grid box
218575850c Andr*0021 C nchp     - integer per-processor number of tiles in tile space
dbae14396f Andr*0022 C lons     - longitude in degrees [im,jm,nSx,nSy]
                0023 C lats     - latitude in degrees [im,jm,nSx,nSy]
614dbd872e Andr*0024 C
                0025 C OUTPUT:
                0026 C
d58f306b82 Jean*0027 C surftype - integer array of land surface types [im,jm,maxtyp,nSx,nSy]
                0028 C tilefrac - real array of corresponding land surface type fractions
                0029 C            [im,jm,maxtyp,nSx,nSy]
                0030 C igrd     - integer array in tile space of grid point number for each
                0031 C            tile [nchp,nSx,nSy]
                0032 C ityp     - integer array in tile space of land surface type for each
                0033 C            tile [nchp,nSx,nSy]
                0034 C chfr     - real array in tile space of land surface type fraction for
                0035 C            each tile [nchp,nSx,nSy]
614dbd872e Andr*0036 C
                0037 C NOTES:
                0038 C       Vegetation type as follows:
                0039 C                  1:  BROADLEAF EVERGREEN TREES
                0040 C                  2:  BROADLEAF DECIDUOUS TREES
                0041 C                  3:  NEEDLELEAF TREES
                0042 C                  4:  GROUND COVER
                0043 C                  5:  BROADLEAF SHRUBS
                0044 C                  6:  DWARF TREES (TUNDRA)
                0045 C                  7:  BARE SOIL
d58f306b82 Jean*0046 C                  8:  DESERT
614dbd872e Andr*0047 C                  9:  GLACIER
                0048 C                 10:  DARK DESERT
                0049 C                100:  OCEAN
                0050 C***********************************************************************
d58f306b82 Jean*0051       IMPLICIT NONE
328acc481c Andr*0052 #include "EEPARAMS.h"
614dbd872e Andr*0053
d58f306b82 Jean*0054       INTEGER myThid,im,jm,maxtyp,nchp,nSx,nSy,Nxg,Nyg
                0055       INTEGER nchptot(nSx,nSy), nchpland(nSx,nSy)
                0056       INTEGER surftype(im,jm,maxtyp,nSx,nSy)
                0057       INTEGER igrd(nchp,nSx,nSy),ityp(nchp,nSx,nSy)
                0058       _RL tilefrac(im,jm,maxtyp,nSx,nSy)
a20b61c7ed Andr*0059       _RL lats(im,jm,nSx,nSy), lons(im,jm,nSx,nSy)
d58f306b82 Jean*0060       _RL chfr(nchp,nSx,nSy),chlt(nchp,nSx,nSy),chlon(nchp,nSx,nSy)
                0061
                0062 C-    local variables:
                0063       CHARACTER*40 vegdata
                0064       INTEGER i,j,k,bi,bj
614dbd872e Andr*0065
d58f306b82 Jean*0066       INTEGER imdata,jmdata,Nxgdata,Nygdata
                0067       INTEGER biglobal,bjglobal
                0068       INTEGER ierr1,kveg
                0069       INTEGER*4 im_32, jm_32, Nxg_32, Nyg_32
                0070       INTEGER*4 iveg_32(im,jm,maxtyp,Nxg,Nyg)
                0071       REAL*4    veg_32(im,jm,maxtyp,Nxg,Nyg)
614dbd872e Andr*0072
d58f306b82 Jean*0073       WRITE(standardMessageUnit,'(2A)') ' FIZHI_INIT_VEG: ',
                0074      &   'defining surface type and fraction: ----------------------'
614dbd872e Andr*0075
26477589cb Jean*0076 #ifdef _BYTESWAPIO
                0077 C     using MDS_BYTESSWAP for sequential acces file does not work
                0078       STOP 'ABNORMAL END: S/R FIZHI_INIT_VEG'
                0079 #endif
d58f306b82 Jean*0080       CALL MDSFINDUNIT( kveg, myThid )
45815db68b Andr*0081       close(kveg)
d58f306b82 Jean*0082       open(kveg,file=vegdata,form='unformatted',access='sequential',
                0083      &                      iostat=ierr1)
614dbd872e Andr*0084       if( ierr1.eq.0 ) then
b71887005c Andr*0085           rewind(kveg)
d58f306b82 Jean*0086           read(kveg)im_32,jm_32,Nxg_32,Nyg_32,iveg_32,veg_32
614dbd872e Andr*0087       else
                0088        print *
                0089        print *, 'Veg Dataset: ',vegdata,' not found!'
                0090        print *
45815db68b Andr*0091        call exit(101)
614dbd872e Andr*0092       endif
45815db68b Andr*0093       close(kveg)
46a396bdaa Andr*0094 #if defined( _BYTESWAPIO ) && defined( ALLOW_MDSIO )
d58f306b82 Jean*0095       CALL MDS_BYTESWAPI4(1,im_32)
                0096       CALL MDS_BYTESWAPI4(1,jm_32)
                0097       CALL MDS_BYTESWAPI4(1,nxg_32)
                0098       CALL MDS_BYTESWAPI4(1,nyg_32)
46a396bdaa Andr*0099 #endif
f9c065a827 Ed H*0100
d58f306b82 Jean*0101       IF (myThid.EQ.1) THEN
                0102        imdata = im_32
                0103        jmdata = jm_32
                0104        Nxgdata = Nxg_32
                0105        Nygdata = Nyg_32
                0106        if( (imdata.ne.im) .or. (jmdata.ne.jm) .or.
                0107      &     (Nxgdata.ne.Nxg) .or. (Nygdata.ne.Nyg) ) then
                0108         print *
                0109         print *, 'Veg Data Resolution is Incorrect! '
                0110         print *,' Model Res: ',im,'x',jm,' Data Res: ',imdata,'x',jmdata
                0111         print *,' Model Nxg Nyg: ',Nxg,' ',Nyg,
                0112      &           ' Data Nxg Nyg: ',Nxgdata,' ',Nygdata
                0113         print *
                0114         call exit(102)
                0115        endif
9a6b9d7b6d Andr*0116       ENDIF
                0117
d58f306b82 Jean*0118       DO bj = myByLo(myThid), myByHi(myThid)
                0119        DO bi = myBxLo(myThid), myBxHi(myThid)
328acc481c Andr*0120
d58f306b82 Jean*0121         biglobal=bi+(myXGlobalLo-1)/im
                0122         bjglobal=bj+(myYGlobalLo-1)/jm
f9c065a827 Ed H*0123 #if defined( _BYTESWAPIO ) && defined( ALLOW_MDSIO )
d58f306b82 Jean*0124         CALL MDS_BYTESWAPR4( im*jm*maxtyp,
                0125      &                       veg_32(1,1,1,biglobal,bjglobal) )
                0126         CALL MDS_BYTESWAPI4( im*jm*maxtyp,
                0127      &                       iveg_32(1,1,1,biglobal,bjglobal) )
be3307cdae Andr*0128 #endif
d58f306b82 Jean*0129         do k = 1,maxtyp
                0130          do j = 1,jm
                0131           do i = 1,im
                0132             surftype(i,j,k,bi,bj) = iveg_32(i,j,k,biglobal,bjglobal)
                0133             tilefrac(i,j,k,bi,bj) = veg_32(i,j,k,biglobal,bjglobal)
                0134           enddo
                0135          enddo
                0136         enddo
614dbd872e Andr*0137
d58f306b82 Jean*0138        ENDDO
f9c065a827 Ed H*0139       ENDDO
                0140
d58f306b82 Jean*0141 C     create chip arrays for :
                0142 C      igrd :  grid index
                0143 C      ityp :  veg. type
                0144 C      chfr :  vegetation fraction
                0145 C      chlon:  chip longitude
                0146 C      chlt :  chip latitude
                0147
                0148 C     nchpland<=nchptot is the actual number of land chips
                0149       WRITE(standardMessageUnit,'(2A)') ' FIZHI_INIT_VEG: ',
                0150      &   'setting surface Tiles:'
                0151
                0152       DO bj = myByLo(myThid), myByHi(myThid)
                0153        DO bi = myBxLo(myThid), myBxHi(myThid)
                0154
                0155 C-       initialise grid index array:
                0156          do i=1,nchp
                0157            igrd(i,bi,bj) = 1
                0158          enddo
                0159
                0160 C-       land points:
                0161          nchpland(bi,bj) = 0
                0162          do k=1,maxtyp
                0163           do j=1,jm
                0164            do i=1,im
                0165              if(surftype(i,j,k,bi,bj).lt.100 .and.
                0166      &            tilefrac(i,j,k,bi,bj).gt.0.) then
                0167                nchpland(bi,bj)  = nchpland(bi,bj) + 1
                0168                igrd (nchpland(bi,bj),bi,bj) = i + (j-1)*im
                0169                ityp (nchpland(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
                0170                chfr (nchpland(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
                0171                chlon(nchpland(bi,bj),bi,bj) = lons(i,j,bi,bj)
                0172                chlt (nchpland(bi,bj),bi,bj) = lats(i,j,bi,bj)
                0173              endif
                0174            enddo
f9c065a827 Ed H*0175           enddo
d58f306b82 Jean*0176          enddo
                0177
                0178 C-       ocean points:
                0179          nchptot(bi,bj) = nchpland(bi,bj)
                0180          do k=1,maxtyp
                0181           do j=1,jm
                0182            do i=1,im
                0183              if(surftype(i,j,k,bi,bj).ge.100 .and.
                0184      &            tilefrac(i,j,k,bi,bj).gt.0.) then
                0185                nchptot(bi,bj)  = nchptot(bi,bj) + 1
                0186                igrd (nchptot(bi,bj),bi,bj) = i + (j-1)*im
                0187                ityp (nchptot(bi,bj),bi,bj) = surftype(i,j,k,bi,bj)
                0188                chfr (nchptot(bi,bj),bi,bj) = tilefrac(i,j,k,bi,bj)
                0189                chlon(nchptot(bi,bj),bi,bj) = lons(i,j,bi,bj)
                0190                chlt (nchptot(bi,bj),bi,bj) = lats(i,j,bi,bj)
                0191              endif
                0192            enddo
f9c065a827 Ed H*0193           enddo
d58f306b82 Jean*0194          enddo
14f15aa407 Andr*0195
d58f306b82 Jean*0196          WRITE(standardMessageUnit,'(2(A,I4),2(A,I10))') '  bi=', bi,
                0197      &    ', bj=', bj, ', # of Land Tiles=', nchpland(bi,bj),
                0198      &                 ', Total # of Tiles=', nchptot(bi,bj)
                0199
                0200        ENDDO
45815db68b Andr*0201       ENDDO
f144613748 Andr*0202
d58f306b82 Jean*0203       WRITE(standardMessageUnit,'(2A)') ' FIZHI_INIT_VEG: done'
                0204
614dbd872e Andr*0205       RETURN
                0206       END