clm_ini.f90

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

      do m = 1, msub
         do l = 1, lpt
            kvec(l,m) = l
            wsg2g(l,m) = 0.
         end do
      end do

      l = 0
      k = 0
      do j = 1, lat
         do i = 1, lon
            if (oroxy(i,j) > 0.) then        !land point
               l = l+1                       !land index
               ixy(l) = i                    !longitude index
               jxy(l) = j                    !latitude index
               ktile(l) = mtile(i,j)
               do m = 1, ktile(l)            !subgrid points
                  k = k + 1                  !subgrid index
                  kvec(l,m) = k              !subgrid index for land point
                  klnd(k) = l                !land index for subgrid point

                  dlat(k) = latixy(i,j)      !latitude in radians
                  dlon(k) = longxy(i,j)      !longitude in radians

                  ivt(k) = surf2d(i,j,m)     !land cover type
                  isc(k) = soic2d(i,j)       !soil color index
                  sand(k) = sand2d(i,j)      !percent of sand
                  clay(k) = clay2d(i,j)      !percent of clay
                  wsg2g(l,m) = wt(i,j,m)     !subgrid weights
               end do
            end if
         end do
      end do

!->   do k = 1, kpt
!->      wetlnd(k) = 0
!->      icetyp(k) = 0
!->      watertyp(k) = 0
!->      if(ivt(k) .eq. 11) wetlnd(k) = 1
!->      if(ivt(k) .eq. 15) icetyp(k) = 1
!->      if(ivt(k) .eq. 17) watertyp(k) = 1
!->   enddo

! ----------------------------------------------------------------------
! initialize time invariant variables, as subgrid vectors of length 
! [kpt], for [clmtc_module]
!                 !please take care the expression, for expample,
!                 !z(2:5,1) is a 1D array with 4 elements;
!                 !z(2:5,1:1) is a 4x1 2D array
! ----------------------------------------------------------------------

      call clmtci (kpt    ,msl    ,ivt    ,isc    ,sand   ,clay   ,&
                   csol   ,porsl  ,phi0   ,bsw    ,dkmg   ,dksatu ,&
                   dkdry  ,albsol ,hksati ,z0m    ,displa ,sqrtdi ,&
                   effcon ,vmax25 ,slti   ,hlti   ,shti   ,hhti   ,&       
                   trda   ,trdm   ,trop   ,gradm  ,binter ,extkn  ,&
                   albvgs ,albvgl ,chil   ,ref    ,tran   ,rootfr ,&
                   z(1,1) ,dz(1,1),zi(0,1))

      write (6,*)
      write (6,*) ('successful set up of clmtc module')
! ----------------------------------------------------------------------
! begin clm initialization
! if continuation run: end of clm initialization because time-varying 
! data in [clmtv_module] has been read in from atmospheric restart 
! file by atmospheric model. atmospheric model also has surface temperature,
! albedos, and snow from restart file.
! ----------------------------------------------------------------------

      dtime = dt_step

      if (nsrest == 0) then            ! initial run

          iyear = mbyear
          jday = mbday
          mcsec = mbsec 

          if(.not. greenwich)then

!         please pay attention when running the off-line regions cases
!         for off-line (stand alone or regions), here the local time 
!         is set as the same as the first grid point

             mcsec = mcsec - longxy(1,1)/(15.*pi/180.)*3600.
             if(mcsec < 0)then
                jday = jday -1
                mcsec = 86400 + mcsec
             endif

             if(jday < 1)then
                iyear = iyear - 1
                if((mod(iyear,4) == 0 .AND. mod(iyear,100) /=0) .OR. mod(iyear,400) == 0)then
                    jday = 366
                 else
                    jday = 365
                 endif
             endif
          endif

!         initialize time-varying data in [clmtv] common block 
!         o rdlsf = false: arbitrary initialization
!         o rdlsf = true : read initial data set

          call clmtvi (kpt   ,msl  ,maxsnl ,lun_ini ,&
                       ivt   ,z0m  ,albsol ,albvgs ,albvgl ,chil  ,ref ,tran, porsl ,&
                       glai  ,gsai ,dlon   ,dlat   ,rdlsf  ,rdlai ,&
                       mcsec ,jday ,xerr   ,zerr ,&
                       dz    ,z    ,zi     ,tss    ,wliq   ,wice ,&
                       rootr ,tlsun,tlsha  ,ldew   ,sag    ,scv  ,snowdp ,&
                       etrc  ,tg   ,albg   ,albv   ,alb    ,ssun ,ssha ,tranc,thermk,extkb,extkd, & 
                       cosz  ,green,fveg   ,fsno   ,sigf   ,lai  ,sai ,&
                       snl)


      else                            ! re-start run
!          read (unit=lun_rst,*) &
!                       mcsec ,jday  ,iyear ,xerr ,zerr ,&
!                       dz    ,z     ,zi    ,tss  ,wliq ,wice ,&
!                       rootr ,tlsun ,tlsha ,ldew ,sag  ,scv  ,snowdp,&
!                       etrc  ,tg    ,albg  ,albv ,alb  ,ssun ,ssha  ,tranc ,thermk ,extkb ,extkd ,& 
!                       cosz  ,green ,fveg  ,fsno ,sigf ,lai  ,sai   ,snl
!
!                       write(6,*)  'balance error of last run =', xerr, zerr
!                       xerr = 0.; zerr = 0.
!                       write(6,*)  mcsec ,jday ,iyear,snl

                        open(888,file='compare.txt',form='formatted')


                        read(unit=lun_rst,*) mcsec, jday, iyear 
                        write(888,*)         mcsec, jday, iyear 
                        read(unit=lun_rst,*) xerr,zerr
                        write(888,*)         xerr,zerr
                        read(unit=lun_rst,*) (( dz(m,n),m=-4,10),n=1,kpt)
                        write(888,*)         dz
                        read(unit=lun_rst,*) ((  z(m,n),m=-4,10),n=1,kpt)
                        write(888,*)         z 
                        read(unit=lun_rst,*) (( zi(m,n),m=-5,10),n=1,kpt) 
                        write(888,*)         zi
                        read(unit=lun_rst,*) ((tss(m,n),m=-4,10),n=1,kpt)
                        write(888,*)         tss
                        read(unit=lun_rst,*) ((wliq(m,n),m=-4,10),n=1,kpt)
                        write(888,*)         wliq
                        read(unit=lun_rst,*) ((wice(m,n),m=-4,10),n=1,kpt)
                        write(888,*)         wliq
                        read(unit=lun_rst,*) (rootr(m),m=1,10)
                        write(888,*)         rootr
                        read(unit=lun_rst,*) tlsun,tlsha
                        write(888,*)         tlsun,tlsha
                        read(unit=lun_rst,*) ldew
                        write(888,*)         ldew
                        read(unit=lun_rst,*) sag,scv,snowdp
                        write(888,*)         sag,scv,snowdp
                        read(unit=lun_rst,*) etrc
                        write(888,*)         etrc
                        read(unit=lun_rst,*) tg
                        write(888,*)         tg
                        read(unit=lun_rst,*) (albg(m),m=1,4)
                        write(888,*)         albg
                        read(unit=lun_rst,*) (albv(m),m=1,4) 
                        write(888,*)         albv
                        read(unit=lun_rst,*) ( alb(m),m=1,4)
                        write(888,*)         alb
                        read(unit=lun_rst,*) (ssun(m),m=1,4) 
                        write(888,*)          ssun 
                        read(unit=lun_rst,*) (ssha(m),m=1,4) 
                        write(888,*)         ssha
                        read(unit=lun_rst,*) thermk 
                        write(888,*)         thermk
                        read(unit=lun_rst,*) extkb,extkd 
                        write(888,*)         extkb,extkd
                        read(unit=lun_rst,*) cosz,green
                        write(888,*)         cosz,green
                        read(unit=lun_rst,*) fveg,fsno,sigf
                        write(888,*)         fveg,fsno,sigf
                        read(unit=lun_rst,*) lai,sai
                        write(888,*)         lai,sai
                        read(unit=lun_rst,*) snl
                        write(888,*)         snl 

                        close(888) 
! for spin-up use
!                      iyear = 1965                          ! for Valdai
!                      iyear = 1986                          ! for Cabauw
!                      iyear = 1986                          ! for Hapex case
!                      iyear = 1983; jday=244; mcsec=14296   ! for arme case
!                      iyear = 1992; jday=1; mcsec=0         ! for abracos case
!                      iyear = 1987; jday=121; mcsec=0       ! for fife case
!                      iyear = 1993; jday=152; mcsec=26659   ! for tucson case 
!                      iyear = 1994                          ! for boreas
      close(lun_rst) 
      endif

      write (6,*)
      write (6,*) ('successful set up of clmtv module')
      print*,'successful set up of clmtv module'
! ----------------------------------------------------------------------
! end clm initialization and return clmtv.h to clmtv_module
! ----------------------------------------------------------------------

! average subgrid albedos, srf temperature, and snow for atmospheric model

! initialize to zero for land only

      do j = 1, lat
         do i = 1, lon
            if (oroxy(i,j) > 0.) then 
               asdirxy(i,j) = 0.
               asdifxy(i,j) = 0.
               aldirxy(i,j) = 0.
               aldifxy(i,j) = 0.
               scv2xy(i,j)  = 0.
               tsxy(i,j)    = 0.
               lwupxy(i,j)  = 0.
            else
               asdirxy(i,j) = -999.
               asdifxy(i,j) = -999.
               aldirxy(i,j) = -999.
               aldifxy(i,j) = -999.
               scv2xy(i,j)  = -999.
               tsxy(i,j)    = -999.
               lwupxy(i,j)  = -999.
            endif
         end do
      end do

! [kpt] vector of subgrid points -> [lpt] vector of land points -> 
! [lon] x [lat] grid 

      do l = 1, lpt         !land point index for [lon] x [lat] grid
         i = ixy(l)         !longitude index
         j = jxy(l)         !latitude index
         do m = 1, ktile(l)
            k = kvec(l,m)   !clm subgrid vector index
            asdirxy(i,j) = asdirxy(i,j) + alb(1,1,k) * wsg2g(l,m)
            asdifxy(i,j) = asdifxy(i,j) + alb(1,2,k) * wsg2g(l,m)
            aldirxy(i,j) = aldirxy(i,j) + alb(2,1,k) * wsg2g(l,m)
            aldifxy(i,j) = aldifxy(i,j) + alb(2,2,k) * wsg2g(l,m)
            scv2xy(i,j)  = scv2xy(i,j)  + scv(k)     * wsg2g(l,m)
            tsxy(i,j)    = tsxy(i,j)    + tg(k) * wsg2g(l,m)
            lwupxy(i,j)  = lwupxy(i,j)  + stefnc*(tg(k)**4) * wsg2g(l,m)
         enddo
      end do

   END SUBROUTINE clm_ini