#include <misc.h>#include <params.h>!----------------------------------------------------------------------- !! BOP!! !MODULE: ice_data!! !DESCRIPTION:	Module to handle dealing with the ICE data.!! Public interfaces:!!	iceini -- Initialization and reading of dataset.!	iceint -- Interpolate dataset ICE to current time.!!----------------------------------------------------------------------- module ice_data!! USES:!  use precision, only: r8  use ice_constants  use pmgrid,    only: plon, plat, masterproc  use ppgrid,    only: pcols, begchunk, endchunk  use phys_grid, only: scatter_field_to_chunk, get_ncols_p,get_lat_all_p  use comsrf,    only: plevmx, icefrac,previcefrac,update_srf_fractions  use physconst, only: tmelt  use commap,    only: clat, clon  implicit none!----------------------------------------------------------------------- ! PUBLIC: Make default data and interfaces private!----------------------------------------------------------------------- !! ! PUBLIC MEMBER FUNCTIONS:!  public iceini   ! Initialization  public iceint   ! Time interpolation of ICE data  logical (kind=log_kind), parameter :: snowice_climatology = .true.!===============================================================================!EOP!===============================================================================!----------------------------------------------------------------------- ! PRIVATE: Everthing else is private to this module!-----------------------------------------------------------------------   private   ! By default all data is private to this module  integer, parameter :: toticesz=2000  real(r8), parameter :: daysperyear = 365.0  ! Number of days in a year  real(r8), allocatable, dimension(:,:,:) :: &      icebdy         ! ICE values on boundary dataset (pcols,begchunk:endchunk,2)  real(r8), allocatable, dimension(:,:) :: &      ice            ! Interpolated model ice values (pcols,begchunk:endchunk)  real(r8) :: cdayicem   ! Calendar day for prv. month ICE values read in  real(r8) :: cdayicep   ! Calendar day for nxt. month ICE values read in        integer :: nm,np   ! Array indices for prv., nxt month ice data  integer :: nm1  integer :: nmshift,npshift ! Array indices for prv., nxt month ice data  integer :: iceid   ! netcdf id for ice variable  integer :: lonsiz  ! size of longitude dimension on ice dataset  integer :: levsiz  ! size of level dimension on ice dataset  integer :: latsiz  ! size of latitude dimension on ice dataset  integer :: timesiz ! size of time dimension on ice dataset  integer :: np1     ! current forward time index of ice dataset  integer :: date_ice(toticesz)! Date on ice dataset (YYYYMMDD)  integer :: sec_ice(toticesz) ! seconds of date on ice dataset (0-86399)   real(r8):: snwbdynh(2)   ! snow height boundary values nrthrn hmsphr  real(r8):: snwbdysh(2)   ! snow height boundary values sthrn hmsphr  real(r8):: snwcnh(12)   ! mean snow cover (m) first of month (nrthrn hmsphr)  real(r8):: snwcsh(12)   ! mean snow cover (m) first of month (sthrn hmsphr)  data snwcnh   /  .23,  .25,  .27,  .29,  .33,  .18, &                 0.,   0.,  .02,  .12,  .18,  .21 /   data snwcsh   /  0.,   0.,  .02,  .12,  .18,  .21, &                  .23,  .25,  .27,  .29,  .33,  .18/!===============================================================================CONTAINS!===============================================================================!======================================================================! PUBLIC ROUTINES: Following routines are publically accessable!======================================================================!----------------------------------------------------------------------- ! ! BOP!! !IROUTINE: iceini!! !DESCRIPTION:!! Initialize the procedure for specifying sea surface temperatures! Do initial read of time-varying ice boundary dataset, reading two! consecutive months on either side of the current model date.!! Method: ! ! Author: L.Bath! !-----------------------------------------------------------------------!! !INTERFACE!subroutine iceini!! !USES:!  use rgrid, only: nlon  use error_messages, only: alloc_err, handle_ncerr  use time_manager, only: get_curr_date, get_curr_calday  use ice_constants#if ( defined SPMD )  use mpishorthand, only: mpicom, mpiint, mpir8#endif!! EOP!!---------------------------Common blocks-------------------------------#include <comctl.h>#include <comlun.h>!---------------------------Local variables-----------------------------  integer dateid                ! netcdf id for date variable  integer secid                 ! netcdf id for seconds variable  integer londimid              ! netcdf id for longitude variable  integer latdimid              ! netcdf id for latitude variable  integer lonid                 ! netcdf id for longitude variable  integer latid                 ! netcdf id for latitude variable  integer timeid                ! netcdf id for time variable  integer nlonid                ! netcdf id for nlon variable (rgrid)  integer cnt3(3)               ! array of counts for each dimension  integer strt3(3)              ! array of starting indices  integer n                     ! indices  integer nlon_ice(plat)        ! number of lons per lat on bdy dataset  integer i                     ! index into chunk  integer j                     ! latitude index  integer k  integer ncol  integer istat                 ! error return  integer lchnk           ! chunk to process  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               ! calendar day (includes yr if no cycling)  real(r8) caldayloc            ! calendar day (includes yr if no cycling)  real(r8) xvar(plon,plat,2)    ! work space !-----------------------------------------------------------------------!! For aqua_planet there is no ice anywhere!  if(aqua_planet)then     icefrac(:pcols,begchunk:endchunk) = 0.0     call update_srf_fractions     return  end if! initialize ice constants	  call init_constants!! Initialize time indices!  nm = 1  np = 2!! Allocate space for data.!  allocate( ice(pcols,begchunk:endchunk), stat=istat )  call alloc_err( istat, 'iceini', 'ice', &       pcols*(endchunk-begchunk+1) )  allocate( icebdy(pcols,begchunk:endchunk,2), stat=istat )  call alloc_err( istat, 'iceini', 'icebdy', &       pcols*(endchunk-begchunk+1)*2 )!! SPMD: Master does all the work.!  if (masterproc) then!! Use year information only if not cycling ice 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     ncdate = yr*10000 + mon*100 + day!! Get and check dimension info!     call wrap_inq_dimid( ncid_sst, 'lon', londimid   )     call wrap_inq_dimid( ncid_sst, 'time', timeid  )     call wrap_inq_dimid( ncid_sst, 'lat', latdimid   )     call wrap_inq_dimlen( ncid_sst, londimid, lonsiz   )     if (lonsiz /= plon) then        write(6,*)'ICEINI: lonsiz=',lonsiz,' must = plon=',plon        call endrun     end if     call wrap_inq_dimlen( ncid_sst, latdimid, latsiz   )     if (latsiz /= plat) then        write(6,*)'ICEINI: latsiz=',latsiz,' must = plat=',plat        call endrun     end if     call wrap_inq_dimlen( ncid_sst, timeid, timesiz   )!! Check to make sure space allocated for time variables is sufficient!     if (timesiz>toticesz) then        write(6,*)'ICEINI:  Allocated space for ice data is insufficient.'        write(6,*)'Please increase parameter toticesz to',timesiz,' and recompile.'        call endrun     end if!! Check to ensure reduced or not grid of dataset matches that of model!     if (fullgrid) then        call wrap_inq_varid( ncid_sst, 'lon', lonid   )     else        call wrap_inq_varid (ncid_sst, 'nlon', nlonid)        call wrap_get_var_int (ncid_sst, nlonid, nlon_ice)        do j=1,plat           if (nlon_ice(j) /= nlon(j)) then              write(6,*)'ICEINI: model grid does not match dataset grid'              call endrun           end if        end do     end if     call wrap_inq_varid( ncid_sst, 'date', dateid   )     call wrap_inq_varid( ncid_sst, 'datesec', secid   )     call wrap_inq_varid( ncid_sst, 'ice_cov', iceid   )     call wrap_inq_varid( ncid_sst, 'lat', latid   )!! Retrieve entire date and sec variables.!     call wrap_get_var_int (ncid_sst,dateid,date_ice)     call wrap_get_var_int (ncid_sst,secid,sec_ice)     if (icecyc) then        if (timesiz<12) then            write(6,*)'ICEINI: ERROR'            write(6,*)'When cycling ice, ice data set must have 12'            write(6,*)'consecutive months of data starting with Jan'           write(6,*)'Current dataset has only ',timesiz,' months'           call endrun        end if        do n = 1,12           if (mod(date_ice(n),10000)/100/=n) then              write(6,*)'ICEINI: ERROR'               write(6,*)'When cycling ice, ice data set must have 12'               write(6,*)'consecutive months of data starting with Jan'              write(6,*)'Month ',n,' of ice data set is out of order'              call endrun           end if        end do     end if     strt3(1) = 1     strt3(2) = 1     strt3(3) = 1     cnt3(1)  = lonsiz     cnt3(2)  = latsiz     cnt3(3)  = 1!! Special code for interpolation between December and January!     if (icecyc) then        n = 12        np1 = 1        call bnddyi(date_ice(n  ), sec_ice(n  ), cdayicem)        call bnddyi(date_ice(np1), sec_ice(np1), cdayicep)        if (caldayloc<=cdayicep .or. caldayloc>cdayicem) then           strt3(3) = n           call wrap_get_vara_realx (ncid_sst,iceid,strt3,cnt3,xvar(1,1,nm))           snwbdynh(nm)=snwcnh(n)           snwbdysh(nm)=snwcsh(n)           strt3(3) = np1                                                 call wrap_get_vara_realx (ncid_sst,iceid,strt3,cnt3,xvar(1,1,np))           snwbdynh(np)=snwcnh(np1)           snwbdysh(np)=snwcsh(np1)           goto 10        end if     end if!! Normal interpolation between consecutive time slices.!     do n=1,timesiz-1        np1 = n + 1        call bnddyi(date_ice(n  ), sec_ice(n  ), cdayicem)        call bnddyi(date_ice(np1), sec_ice(np1), cdayicep)        if (.not.icecyc) then           yr = date_ice(n)/10000           cdayicem = cdayicem + yr*daysperyear           yr = date_ice(np1)/10000           cdayicep = cdayicep + yr*daysperyear        end if        if (caldayloc>cdayicem .and. caldayloc<=cdayicep) then           strt3(3) = n           call wrap_get_vara_realx (ncid_sst,iceid,strt3,cnt3,xvar(1,1,nm))           snwbdynh(nm)=snwcnh(n)           snwbdysh(nm)=snwcsh(n)           strt3(3) = np1                                                 call wrap_get_vara_realx (ncid_sst,iceid,strt3,cnt3,xvar(1,1,np))           snwbdynh(np)=snwcnh(np1)           snwbdysh(np)=snwcsh(np1)           goto 10        end if     end do     write(6,*)'ICEINI: Failed to find dates bracketing ncdate, ncsec=', ncdate, ncsec     call endrun10   continue     write(6,*)'ICEINI: Read ice data for dates ',date_ice(n),sec_ice(n), &          ' and ',date_ice(np1),sec_ice(np1)#if (defined SPMD )     call mpibcast( timesiz, 1, mpiint, 0, mpicom )     call mpibcast( date_ice, toticesz, mpiint, 0, mpicom )     call mpibcast( sec_ice, toticesz, mpiint, 0, mpicom )     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 )  else     call mpibcast( timesiz, 1, mpiint, 0, mpicom )     call mpibcast( date_ice, toticesz, mpiint, 0, mpicom )     call mpibcast( sec_ice, toticesz, mpiint, 0, mpicom )