ice_data.f90

来自「CCSM Research Tools: Community Atmospher」· F90 代码 · 共 664 行 · 第 1/2 页

     call mpibcast( cdayicem, 1, mpir8, 0, mpicom )
     call mpibcast( cdayicep, 1, mpir8, 0, mpicom )
     call mpibcast( np1, 1, mpiint, 0, mpicom )
     call mpibcast( snwbdynh, 2, mpir8, 0, mpicom )
     call mpibcast( snwbdysh, 2, mpir8, 0, mpicom )
#endif
  end if

  call scatter_field_to_chunk(1,1,2,plon,xvar,icebdy)
!
! Uncomment the following to get fractional model to give bit for bit 
! with non fractional
!
!!$  do lchnk=begchunk,endchunk
!!$     ncol = get_ncols_p(lchnk)
!!$     do i=1,ncol
!!$        do k=1,2
!!$           if (icebdy(i,lchnk,k)<=tsice) then
!!$              icebdy(i,lchnk,k) = 1.0
!!$           else
!!$              icebdy(i,lchnk,k) = 0.
!!$           end if
!!$        end do
!!$     end do
!!$  end do

  return
end subroutine iceini

!----------------------------------------------------------------------- 
! 
! BOP
!
! !IROUTINE: iceint
!
! !DESCRIPTION:
!
! if "aqua_planet", specify ICE's analytically (Jerry Olson).
! Otherwise, time interpolate ICE's to current time, reading in new monthly data if
! necessary.
!
! Method: 
! 
! Author: L.Bath
! 
!-----------------------------------------------------------------------
!
! !INTERFACE:
!
subroutine iceint
!
! !USES:
!
  use rgrid, only: nlon
  use comsrf, only: snowhice, sicthk, icefrac, previcefrac,surface_state2d
  use time_manager, only: is_first_step, get_curr_date, get_curr_calday
  use ice_dh, only: prognostic_icesnow
#if ( defined SPMD )
  use mpishorthand, only: mpicom, mpir8
#endif
!
! EOP
!
!---------------------------Common blocks-------------------------------
#include <comctl.h>
#include <comlun.h>
!---------------------------Local variables-----------------------------
  integer cnt3(3)        ! array of counts for each dimension
  integer strt3(3)       ! array of starting indices
  integer i,j,lchnk      ! indices
  integer ncol           ! number of columns in current chunk
  integer ntmp           ! temporary
  real(r8) fact1, fact2  ! time interpolation factors
  integer :: yr, mon, day! components of a date
  integer :: ncdate      ! current date in integer format [yyyymmdd]
  integer :: ncsec       ! current time of day [seconds]
  real(r8) :: calday     ! current calendar day
  real(r8) caldayloc     ! calendar day (includes yr if no cycling)
  real(r8) deltat        ! time (days) between interpolating ice data
!
! Aqua planet variables
!
  real(r8) pi            ! 3.14159...
  real(r8) pio180        ! pi/180.
  real(r8) tmp           ! temporary
  real(r8) tmp1          ! temporary
  real(r8) t0_max        ! max reference temperature
  real(r8) t0_min        ! min reference temperature
  real(r8) t0_max6       ! max asymmetric reference temperature for option 6
  real(r8) t0_max7       ! max asymmetric reference temperature for option 7
  real(r8) maxlat        ! cutoff latitude poleward of which ICE = 0 deg C
  real(r8) shift         ! number of degrees peak ICE is shifted off equator
  real(r8) shift9        ! number of degrees peak ICE is shifted off equator for opt. 9
  real(r8) shift10       ! number of degrees peak ICE is shifted off equator for opt. 10
  real(r8) latcen        ! center of asymmetric ICE forcing
  real(r8) latrad6       ! radius of asymmetric ICE forcing for option 6
  real(r8) latrad8       ! radius of asymmetric ICE forcing for option 8
  real(r8) loncen        ! center of asymmetric ICE forcing
  real(r8) lonrad        ! radius of asymmetric ICE forcing
  real(r8) xvar(plon,plat,2)    ! work space 
  integer  ice_option    ! option of analytical ICE algorithm
  integer :: lats(pcols)           ! array of latitude indices
!
!-----------------------------------------------------------------------
     if(aqua_planet) return
!
! SPMD: Master does all the work.  Sends needed info to slaves
!
!
! Use year information only if a multiyear dataset
!
     calday = get_curr_calday()
     call get_curr_date(yr, mon, day, ncsec)
     if (icecyc) then
        caldayloc = calday
     else
        caldayloc = calday + yr*daysperyear
     end if

     if (masterproc) then

        strt3(1) = 1
        strt3(2) = 1
        strt3(3) = 1
        cnt3(1)  = lonsiz
        cnt3(2)  = latsiz
        cnt3(3)  = 1

     endif
!
! If model time is past current forward ice timeslice, read in the next
! timeslice for time interpolation.  Messy logic is for icecyc = .true. 
! interpolation between December and January (np1==1).  Note that 
! np1 is never 1 when icecyc is .false.
!
     if (caldayloc > cdayicep .and. .not. (np1==1 .and. caldayloc>cdayicem)) then
        if (icecyc) then
           np1 = mod(np1,12) + 1
        else
           np1 = np1 + 1
        end if
        if (np1>timesiz) then
           if (masterproc) then
              write(6,*)'ICEINT: Attempt to read past end of ICE dataset'
           endif
           call endrun
        end if
        cdayicem = cdayicep
        call bnddyi(date_ice(np1), sec_ice(np1), cdayicep)

        if (.not.icecyc) then
           yr = date_ice(np1)/10000
           cdayicep = cdayicep + yr*daysperyear
        end if

        if (np1==1 .or. caldayloc<=cdayicep) then
           ntmp = nm
           nm = np
           np = ntmp
           if (masterproc) then
              strt3(3) = np1
              call wrap_get_vara_realx (ncid_sst,iceid,strt3,cnt3,xvar(1,1,np))
              write(6,*)'ICEINT: Read ice for date (yyyymmdd) ',date_ice(np1), &
                   ' sec ',sec_ice(np1)
           endif
           snwbdynh(np)=snwcnh(np1)
           snwbdysh(np)=snwcsh(np1)
           call scatter_field_to_chunk(1,1,1,plon,xvar(1,1,np),icebdy(1,begchunk,np))
        else
           if (masterproc) then
              write(6,*)'ICEINT: Input ice for date',date_ice(np1), &
                   ' sec ',sec_ice(np1), 'does not exceed model date',ncdate,&
                   ' sec ',ncsec,' Stopping.'
           endif
           call endrun
        end if
     end if
!
! Time interpolation.  Account for December-January interpolation if
! cycling ice dataset.  Again note that np1 is never 1 when icecyc is false
!
     if (np1==1) then                    ! Dec-Jan interpolation
        deltat = cdayicep + daysperyear - cdayicem
        if (caldayloc>cdayicep) then      ! We're in December
           fact1 = (cdayicep + daysperyear - caldayloc)/deltat
           fact2 = (caldayloc - cdayicem)/deltat
        else                                ! We're in January
           fact1 = (cdayicep - caldayloc)/deltat
           fact2 = (caldayloc + daysperyear - cdayicem)/deltat
        end if
     else
        deltat = cdayicep - cdayicem
        fact1 = (cdayicep - caldayloc)/deltat
        fact2 = (caldayloc - cdayicem)/deltat
     end if
!
! Check sanity of time interpolation calculation to within 32-bit roundoff
!
     if (abs(fact1+fact2-1.)>1.e-6 .or. &
             fact1>1.000001 .or. fact1<-1.e-6 .or. &
             fact2>1.000001 .or. fact2<-1.e-6) then 
        if (masterproc) then
           write(6,*)'ICEINT: Bad fact1 and/or fact2=',fact1,fact2
        endif
        call endrun
     end if

     previcefrac(:,:) = icefrac(:,:)
     do lchnk=begchunk,endchunk
        ncol = get_ncols_p(lchnk)
        call get_lat_all_p(lchnk, ncol, lats)
        do i=1,ncol
              icefrac(i,lchnk) = icebdy(i,lchnk,nm)*fact1 + icebdy(i,lchnk,np)*fact2
              icefrac(i,lchnk)=min(icefrac(i,lchnk),1.)
              icefrac(i,lchnk)=max(icefrac(i,lchnk),0.)
!
! If we are not prognosing snow then use a snow climatology.
! Compute snow height for given snow climatology phase shift snow
! height 6 months for southern hemisphere this only works when
! using ice climatology dataset
!
              if (.not. prognostic_icesnow) then
                 if (.not.icecyc) then
                    write(6,*)'ice_data: Snowice climatology option ',&
                         'not valid for icecyc = .false.'
                    call endrun
                 end if
                 if (clat(lats(i)).gt.0) then
                    snowhice(i,lchnk)=snwbdynh(nm)*rhos/rhofresh * fact1 &
                         + snwbdynh(np)*rhos/rhofresh * fact2
                 else 
                    snowhice(i,lchnk)=snwbdysh(nm)*rhos/rhofresh*fact1 + &
                         snwbdysh(np)*rhos/rhofresh*fact2
                 end if
              end if
        end do
     end do
     if (is_first_step()) previcefrac(:,:) = icefrac(:,:)
!
! get new ocean/ice percentages since iceint will change ice coverage
!
     call update_srf_fractions
!
! check for new ice and give it some default values for sicthk and snowhice
! and tssub

     call newiceproperties
     return
   end subroutine iceint

!======================================================================
! PRIVATE ROUTINES: Following routines are Private to this module
!======================================================================
!----------------------------------------------------------------------- 
! 
! BOP
!
! !IROUTINE: newiceproperties
!
! Purpose: 
!	Set the initial sea-ice surface variables associated with ICE.
!	This routine is private to this module.
!
!-----------------------------------------------------------------------
!
! !INTERFACE:
!
subroutine newiceproperties
!
! !USES:

  use comsrf, only: snowhice, sicthk, icefrac, previcefrac,surface_state2d,tsice
  use sst_data, only: sst
  use time_manager, only: is_first_step
  use ice_constants, only: Tfrezk
  use ice_dh, only: reset_csim_iceprops
!
! EOP
!
!---------------------------Local variables-----------------------------
  integer lchnk, i,k  ! loop indices
  integer ncol        ! number of columns in current chunk
  real(r8) tsurf
  integer :: lats(pcols)           ! array of latitude indices
!
! Set initial ice surface values
!
!Set constant sea ice thickness for each hemisphere. For new ice and
!for the reset_csim_ice_props=.T. namelist variable set
!initial temperature profile to freezing

! if reset_csim_iceprops set the previous ice fraction to 0 to make
! all ice appear as though it is new sea ice. This will effectively
! reset all ice properties to a the base state of new ice.  This is
! only provided on the first timestep of a run.

  if (is_first_step().and. reset_csim_iceprops) then
     write(6,*)'**********Resetting previcefrac to 0********'
     previcefrac(:,:)=0.	
  end if

  do lchnk=begchunk,endchunk
     ncol = get_ncols_p(lchnk)
     call get_lat_all_p(lchnk, ncol, lats)
     do i=1,ncol
        if (icefrac(i,lchnk) > 0._r8) then
           if (clat(lats(i)).gt.0) then
              sicthk(i,lchnk) = 2.0
           else
              sicthk(i,lchnk) = 1.0
           endif
           if ( previcefrac(i,lchnk)==0.) then !newice
	      snowhice(i,lchnk) = 0._r8
              tsice(i,lchnk) = Tfrezk
              do k=1,plevmx
                 surface_state2d(lchnk)%tssub(i,k)=Tfrezk
              end do
           end if
        else !tssub for non ice boxes filled with open ocean temp
           snowhice(i,lchnk) = 0._r8
           sicthk(i,lchnk) = 0._r8
           tsice(i,lchnk) = 0._r8
           do k=1,plevmx
              surface_state2d(lchnk)%tssub(i,k)=sst(i,lchnk)+273.16
           end do
        end if
     end do
  end do
  return
end subroutine newiceproperties
end module ice_data