clm_ini.f90

来自「CLM集合卡曼滤波数据同化算法」· F90 代码 · 共 493 行 · 第 1/2 页

   subroutine clm_ini (lon     ,lat     ,msub ,&
                       nsrest  ,rdlsf   ,rdlai ,&
                       dt_step ,mbyear  ,mbday   ,mbsec ,greenwich ,&
                       lun_ini ,lun_rst ,lun_srf ,lun_lai)

! ======================================================================
!      Source file: clm_ini.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! Description:
! initialization routine for land surface model. initialize land surface
! variables for atmosphere model if not a continutation run.
! calling sequence is:

! clm_ini:            clm initialization
!   -> clm2d_allocate:allocates memory for CLM 2D variables 
!   -> clm1d_allocate:allocates memory for CLM 1D variables 
!   -> clmtci:        set time invariant variables for clmtc_module
!   -> clmtvi:        set time varying variables for clmtv_module
!      -> snowini:    snow layer number, thickness and depth initialization
!      -> fractsnow: fraction of snow cover
!      -> clmzen:     solar zenith angle
!      -> albedo:     surface albedo

! the land surface model works by gathering all the land points on a
! [lon] x [lat] grid into a vector of [lpt] land points. this is
! then expanded into a vector of [kpt] subgrid points, allowing for
! up to [msub] subgrid points per land point. the [kpt] subgrid points
! are processed in strips of [numpnt] points for [numlv] calls to the
! land surface model. this subroutine: 

!   o reads time-invariant surface data on [lon] x [lat] grid 
!   o reads initial data  
!   o builds the appropriate subgrid <-> grid mapping indices  
!   o initializes time constant [clmtc_module] and time variant [clmtv_module]
!   o if initial run: provides required land surface variables to atm model
! ======================================================================

   USE phycon_module
   USE enkfpar_module
   use clmpar_module
   USE clmtc_module
   USE clmtv_module
   USE clm2d_module

   IMPLICIT NONE

! ------------------------ input/output variables ----------------------
! input from atmospheric model

   integer, INTENT(in) :: &
         lon,             &!longitude points
         lat,             &!latitude points
         msub,            &!number of submesh
         mbyear,          &!base year of run 
         mbday,           &!base day of run 
         mbsec,           &!base seconds of base date 

         nsrest,          &!0 = initial run, > 0 = continuation run

         lun_srf,         &!logical unit number of surface data
         lun_lai,         &!logical unit number of LAI data
         lun_ini,         &!logical unit number of initial data
         lun_rst           !logical unit number of restart file

   logical, intent(in) :: &
         rdlsf,           &!true if read in initial data set
         rdlai,           &!true if read in LAI data set
         greenwich         !true if starting time is greenwich time
                           !false -> local time

   real, INTENT(in) :: dt_step !time step (seconds)

! output to atmospheric model
!* real, dimension(lon,lat), INTENT(out) :: & 
!*       asdirxy,         &!albedo - visible waveband, direct
!*       asdifxy,         &!albedo - visible waveband, diffuse
!*       aldirxy,         &!albedo - near infrared waveband, direct
!*       aldifxy,         &!albedo - near infrared waveband, diffuse
!*       oroxy,           &!0 = non-land. 1 = land
!*       tsxy,            &!surface (skin) temperature
!*       scv2xy,          &!water equivalent snow (mm)
!*       lwupxy            !longwave up flux (W/m**2)

! ------------------------ local variables -----------------------------
! run control parameters inlcuded as [clmexp] namelist. except as follows,
! these are in the [clmtc] common block

   integer i,j,k,l,m       !indices

   integer &
         surf2d(lon,lat,msub), &!vegetation type
         soic2d(lon,lat),      &!soil color
         mtile(lon,lat)         !number of tiles in grid

   real  sand2d(lon,lat),      &!percent sand
         clay2d(lon,lat),      &!percent clay
         wt (lon,lat,msub),    &!subgrid weights
         latixy(lon,lat),      &!latitude in radians
         longxy(lon,lat)        !longitude in radians

   integer &
         ireg,            &!region with the order of 80-90S, 70-80S, ... 
         iveg              !IGBP vegetation type (1 to 18)

   real  silt2d,          &!percent silt
         pi,              &!pie
         zoroxy            !temporay use 

   integer ::m,n 

! ----------------------------------------------------------------------
! read time-invariant boundary data on [lon] x [lat] grid.

!    o first [lat] values: number of longitude points for each latitude.
!      this allows for variable longitudinal resolution for each latitude

! remaining data is for each grid cell:

!    o 1st col: latitude  at center of grid cell (degrees)
!    o 2nd col: longitude at center of grid cell (degrees)
!    o 3th col: soil color (1 to 8) for use with soil albedos
!    o 4th col: soil texture, %sand, for thermal and hydraulic properties
!    o 5th col: soil texture, %silt, for thermal and hydraulic properties
!    o 6th col: soil texture, %clay, for thermal and hydraulic properties
!    o 7th col: number of tiles in grid (i,j)
!    o 8th col: surface type, for use as multiple subgrid point
!    o 9th col: subgrid weight
!     (tile 1: dominant type 1, tile 2: dominent type 2, tile 3: bare soil,
!      tile 4: perminent wetlands, tile 5: lake); ocean type index: 0 and 
!              no mosaic used for ocean.
!      if number show type index = 0 and fraction = 0, mean: no submesh defined
!
! ------------------ note ------------------
! the model required the same longitudinal resolution for 
! all latitude strip. For the variable longitudinal resolution
! cases, please assign the surface types and soil character:
! 0 to ocean grids and -999 to the land grids 
! which are not included in the calcultion.
! ----------------------------------------------------------------------

! Note: for coupling run the 2D allocations have to be assigned
! before calling the land component, an alternative is that remove the
! clm2d_module and clm2d_allocation, and the assign the dimension in the 
! GCM interface subroutine.

      call  clm2d_allocate (lon, lat, msl)

      do j = 1, lat
         do i = 1, lon
            read (lun_srf,*) latixy(i,j),longxy(i,j),&
             soic2d(i,j),sand2d(i,j),silt2d,clay2d(i,j),&
             mtile(i,j),(surf2d(i,j,k),k=1,msub),(wt(i,j,k),k=1,msub)
         end do
      end do
      close (lun_srf)

! convert latitudes and longitudes from degress to radians
      pi = 4.*atan(1.)                !pi
      do j = 1, lat
         do i = 1, lon
            latixy(i,j) = latixy(i,j)*pi/180. 
            longxy(i,j) = longxy(i,j)*pi/180. 
         end do
      end do

! return 2-d orography based on surf2d
! please pay attention on the coupling run, we assign "zero" to all SEA grid 
! (including "sea-ice"), when you do the re-start run, please make sure the 
! following assignment should not change the sea-ice definition.

      do j = 1, lat
         do i = 1, lon
            zoroxy = 0.
            do k = 1, mtile(i,j)
               zoroxy = zoroxy + float(surf2d(i,j,k))
            enddo
            if (zoroxy > 0.) then
               oroxy(i,j) = 1.
            else
               oroxy(i,j) = 0.
            end if
         end do
      end do

! ----------------------------------------------------------------------
! read leaf and stem area index
! imon: month (1 -> 12)
! ireg: each region with the order of 90-80S, 80-70S, ..., 80 - 90N (1 -> 18)
! iveg: IGBP land cover type (1 -> 18)
! glai: LAI
! gsai: stem and dead leaf area index 
! ----------------------------------------------------------------------
      if(rdlai)then
         do i = 1, 18            !land cover type
            do j = 1, 18         !region
               do k = 1, 12      !month
                  read (lun_lai,*) iveg, ireg, glai(i,j,k), gsai(i,j,k)
               enddo
            enddo
         enddo
         close (lun_lai)
      else
         write(6,*) 
         write(6,*) 'waiting to creat LAI'
      endif

! ----------------------------------------------------------------------
! build mapping indices and weights: [lon] x [lat] grid <->
! [lpt] vector of land points <-> [kpt] vector of subgrid points. 

! example: [kpt] vector of subgrid points -> [lon] x [lat] grid
! for 1 <= l <= lpt: 
! t_xy( ixy(l), jxy(l) ) =  t_vec( kvec(l,   1) ) * wsg2g(l,   1) + ... +
!                           t_vec( kvec(l,msub) ) * wsg2g(l,msub) 

! if the grid cell has less than [msub] subgrid points, an arbitrary 
! value of [kvec] is used to index the subgrid vector. the weight is zero

! example: [lon] x [lat] grid -> [kpt] vector of subgrid points
! for 1 <= k <= kpt: t_vec(k) = t_xy( ixy(klnd(k)), jxy(klnd(k)) )
! --------------------------------------------------------------------
      
      k = 0
      l = 0
      do j = 1, lat
         do i = 1, lon
               if (oroxy(i,j) > 0.) then             !land point
                   l = l+1                           !land index
                   do m = 1, mtile(i,j)
                      if (wt(i,j,m) .gt. 0.) then    !valid subgrid point
                          k = k+1                    !subgrid index
                      end if
                   end do
               end if
         end do
      end do

      lpt = l
      kpt = k

      call clm1d_allocate (lpt, kpt, msub, maxsnl, msl)

! ivt, isc, sand, and clay. 
! initialize [kvec] for [msub] subgrid types for [lpt] land cells
! to arbitrary value with weight = 0