#include <misc.h>#include <params.h>subroutine oznini!----------------------------------------------------------------------- ! ! Purpose: Do initial read of time-variant ozone boundary dataset, containing!          ozone mixing ratios as a function of latitude and pressure.  Read two!          consecutive months between which the current date lies.  Routine!          RADOZ2 then evaluates the two path length integrals (with and without!          pressure weighting) from zero to the interfaces between the input!          levels.  It also stores the contribution to the integral from each!          layer.! ! Method: Call appropriate netcdf wrapper routines and interpolate to model grid! ! Author: CCM Core Group! !-----------------------------------------------------------------------   use precision   use pmgrid   use comozp   use pspect   use rgrid   use commap   use time_manager, only: get_curr_date, get_perp_date, get_curr_calday, &                           is_perpetual#if ( defined SPMD )   use mpishorthand#endif   implicit none#include <comctl.h>#include <comlun.h>!-----------------------------------------------------------------------   include 'netcdf.inc'!-----------------------------------------------------------------------!! Local workspace!   integer dateid                          ! netcdf id for date variable   integer secid                           ! netcdf id for seconds variable   integer londimid                        ! netcdf id for longitude dimension   integer latdimid                        ! netcdf id for latitude dimension   integer levdimid                        ! netcdf id for level dimension   integer lonid                           ! netcdf id for longitude variable   integer latid                           ! netcdf id for latitude variable   integer levid                           ! netcdf id for level variable   integer timeid                          ! netcdf id for time variable   integer dimids(nf_max_var_dims)         ! variable shape   integer cnt4(4)                         ! array of counts for each dimension   integer strt4(4)                        ! array of starting indices   integer i, k, lat, n                    ! longitude, level, latitude, time indices   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), allocatable :: tmpozmix(:,:)  ! temporary ozmix array   real(r8), allocatable :: oznbdym(:,:,:) ! ozone data previous time sample   real(r8), allocatable :: oznbdyp(:,:,:) ! ozone data next time sample   nm = 1   np = 2!! SPMD: Master does all the work.  Sends needed info to slaves!   if (masterproc) then!! Use year information only if not cycling ozone dataset!      calday = get_curr_calday()      if ( is_perpetual() ) then         call get_perp_date(yr, mon, day, ncsec)      else         call get_curr_date(yr, mon, day, ncsec)      end if      ncdate = yr*10000 + mon*100 + day      if (ozncyc) then         caldayloc = calday      else         caldayloc = calday + yr*365.      end if!! Get and check dimension info!      CALL WRAP_INQ_DIMID( ncid_oz, 'lon', londimid   )      CALL WRAP_INQ_DIMID( ncid_oz, 'lev', levdimid   )      CALL WRAP_INQ_DIMID( ncid_oz, 'time', timeid  )      CALL WRAP_INQ_DIMID( ncid_oz, 'lat', latdimid   )      CALL WRAP_INQ_DIMLEN( ncid_oz, londimid, lonsiz   )      CALL WRAP_INQ_DIMLEN( ncid_oz, levdimid, levsiz   )      CALL WRAP_INQ_DIMLEN( ncid_oz, latdimid, latsiz   )      CALL WRAP_INQ_DIMLEN( ncid_oz, timeid, timesiz   )      CALL WRAP_INQ_VARID( ncid_oz, 'date', dateid   )      CALL WRAP_INQ_VARID( ncid_oz, 'datesec', secid   )      CALL WRAP_INQ_VARID( ncid_oz, 'OZONE', oznid   )      CALL WRAP_INQ_VARID( ncid_oz, 'lon', lonid   )      CALL WRAP_INQ_VARID( ncid_oz, 'lat', latid   )      CALL WRAP_INQ_VARID( ncid_oz, 'lev', levid   )      CALL WRAP_INQ_VARDIMID (ncid_oz, oznid, dimids)      if (dimids(1) /= londimid .or. dimids(2) /= levdimid .or. dimids(3) /= latdimid) then         write(6,*)'OZNINI: Data must be ordered lon, lev, lat, time'         call endrun      end if!! Dynamically allocated memory for module comozp !      allocate (ozlon(lonsiz))      allocate (ozlat(latsiz))      allocate (date_oz(timesiz))      allocate (sec_oz(timesiz))      allocate (pin(levsiz))      allocate (ozmixm(plond,levsiz,plat,2))      allocate (ozmix(plond,levsiz,plat))!! Locally dynamic that will be deallocated before "return"!      allocate (oznbdym(lonsiz,levsiz,latsiz))      allocate (oznbdyp(lonsiz,levsiz,latsiz))      allocate (tmpozmix(levsiz,plat))!! Retrieve longitude, latitude and level arrays for interpolation.!      CALL WRAP_GET_VAR_REALX (NCID_OZ, lonid,ozlon)      CALL WRAP_GET_VAR_REALX (NCID_OZ, latid,ozlat)      CALL WRAP_GET_VAR_REALX (NCID_OZ, levid,pin)!! Convert from millibars to pascals!      do k=1,levsiz         pin(k) = pin(k)*100.      end do!! Retrieve entire date and sec variables.!      CALL WRAP_GET_VAR_INT (ncid_oz,dateid,date_oz)      CALL WRAP_GET_VAR_INT (ncid_oz,secid,sec_oz)      if (ozncyc) then         if (timesiz < 12) then             write(6,*)'OZNINI: When cycling ozone, dataset must have 12 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_oz(n),10000)/100 /= n) then               write(6,*)'OZNINI: When cycling ozone, dataset must have 12 consecutive ', &                         'months of data starting with Jan'               write(6,*)'Month ',n,' of dataset says date=',date_oz(n)               call endrun            end if         end do      end if      strt4(1) = 1      strt4(2) = 1      strt4(3) = 1      cnt4(1)  = lonsiz      cnt4(2)  = levsiz      cnt4(3)  = latsiz      cnt4(4)  = 1!! Special code for interpolation between December and January!      if (ozncyc) then         n = 12         np1 = 1         call bnddyi(date_oz(n  ), sec_oz(n  ), cdayozm)         call bnddyi(date_oz(np1), sec_oz(np1), cdayozp)         if (caldayloc <= cdayozp .or. caldayloc > cdayozm) then            strt4(4) = n            call wrap_get_vara_realx (ncid_oz,oznid,strt4,cnt4,oznbdym)            strt4(4) = np1            call wrap_get_vara_realx (ncid_oz,oznid,strt4,cnt4,oznbdyp)            goto 10         end if      end if!! Normal interpolation between consecutive time slices.!      do n=1,timesiz-1         np1 = n + 1         call bnddyi(date_oz(n  ), sec_oz(n  ), cdayozm)         call bnddyi(date_oz(np1), sec_oz(np1), cdayozp)         if (.not.ozncyc) then            yr = date_oz(n)/10000            cdayozm = cdayozm + yr*365.            yr = date_oz(np1)/10000            cdayozp = cdayozp + yr*365.         end if         if (caldayloc > cdayozm .and. caldayloc <= cdayozp) then            strt4(4) = n            call wrap_get_vara_realx (ncid_oz,oznid,strt4,cnt4,oznbdym)            strt4(4) = np1            call wrap_get_vara_realx (ncid_oz,oznid,strt4,cnt4,oznbdyp)            goto 10         end if      end do      write(6,*)'OZNINI: Failed to find dates bracketing ncdate, ncsec=', ncdate, ncsec      call endrun10    continue      write(6,*)'OZNINI: Read ozone data for dates ',date_oz(n), &                sec_oz(n),' and ',date_oz(np1),sec_oz(np1)!! Spatial interpolation.  If ozone dataset is 2-d (i.e. lonsiz = 1) and ! thus only latitude interpolation is necessary, expand to 3-d after ! interpolation.!      if (lonsiz == 1) then         call lininterp (oznbdym ,ozlat   ,levsiz  ,latsiz  ,tmpozmix, &                         latdeg  ,plat    )         do lat=1,plat            do k=1,levsiz               do i=1,nlon(lat)                  ozmixm(i,k,lat,nm) = tmpozmix(k,lat)               end do            end do         end do         call lininterp (oznbdyp ,ozlat   ,levsiz  ,latsiz  ,tmpozmix, &                         latdeg  ,plat)         do lat=1,plat            do k=1,levsiz               do i=1,nlon(lat)                  ozmixm(i,k,lat,np) = tmpozmix(k,lat)               end do            end do         end do      else         call bilin (oznbdym, ozlon, ozlat, lonsiz, lonsiz, &                     levsiz, levsiz, latsiz, ozmixm(1,1,1,nm), londeg, &                     latdeg, plond, nlon, levsiz, plat)         call bilin (oznbdyp, ozlon, ozlat, lonsiz, lonsiz, &                     levsiz, levsiz, latsiz, ozmixm(1,1,1,np), londeg, &                     latdeg, plond, nlon, levsiz, plat)      end if!! Deallocate dynamic memory for local workspace.  NOT for pointers in common.!      deallocate (oznbdym)      deallocate (oznbdyp)      deallocate (tmpozmix)   end if#if ( defined SPMD )   call mpibcast (levsiz, 1, mpiint, 0, mpicom)   if (.not.masterproc) then      allocate (pin(levsiz))      allocate (ozmix(plond,levsiz,plat))   end if   call mpibcast (pin, levsiz, mpir8, 0, mpicom)#endif   returnend subroutine oznini