mixed layer depth ...

Colleagues

Here is code in R to calculate mixed layer depth. The code is part of
the rcalcofi package written by my colleague Ed Weber at NOAA SWFSC.

The method is based on Kara et al (2000) "An optimal definition for
ocean mixed layer depth". JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 105, NO.
C7, PAGES 16,803-16,821.
--------------------------------------------------
`calculate.mld` <-
function(sigma, z, deltaSigma = 0.125)
   {
     # check that the number of sigma and z values are equal
     if (length(sigma) != length(z))
     {
         stop('sigma and z vectors must have equal length')
     }

     # remove the na's
     keepers <- !(is.na(z) | is.na(sigma))
     z <- z[keepers]
     sigma <- sigma[keepers]

     # return an NA and a warning if there aren't at least two
     # numeric values
     if (length(sigma) < 2)
     {
         warning('fewer than two valid sigma and depth-value
combinations entered,
                 NA returned')
         NA
     } else {
         # I use negative depths to be consistent with the Scripps database
         if (all(z >= 0))
         {
             z <- z * -1
             pos <- TRUE
         } else {
             pos <- FALSE
         }

         # be sure the data are sorted by depths
         ord <- order(z, decreasing = TRUE)
         z <- z[ord]
         sigma <- sigma[ord]
         #z <- sort(z, decreasing = TRUE)
         #sigma <- sigma[order(z, decreasing = TRUE)]

         # Manuscript uses a z of 10 m as the initial sigmaRef, but we will
         # take the closest value in case the 10-m measurement is missing.
         minDepth <- which(abs(z + 10) == min(abs(z + 10)))
         minDepth <- ifelse(length(minDepth) > 1, minDepth[2], minDepth)

         sigmaRef <- sigma[minDepth]
         sigma <- sigma[minDepth:length(sigma)]
         z <- z[minDepth:length(z)]
         diffs <- abs(sigma - sigmaRef)

         # if sigma never changes by at least deltaSigma, return the
lowest depth
         # Otherwise proceed
         if (max(diffs, na.rm = TRUE) >= deltaSigma)
         {
             # the uniform region, if present, occurs where the change
between
             # any two points is <= deltaSigma * 1/10, and the change in
sigma over
             # the profile has not yet exceeded deltaSigma
             uniformRegion <- (abs(sigma[2:length(sigma)] -
                                   sigma[1:(length(sigma) - 1)]) <=
                                   (deltaSigma / 10)) &
                              (diffs[2:length(diffs)] < deltaSigma)
             if (any(uniformRegion))
             {
                 sigmaRefPos <- max(which(uniformRegion))
                 # change sigmaRef to the base of the uniform region
                 sigmaRef <- sigma[sigmaRefPos]
                 # calculate change from the new reference
                 reachedDeltaSigma <-
                     abs(sigma[(sigmaRefPos + 1):length(sigma)] -
sigmaRef) >=
                     deltaSigma
                 # if any deeper measurements of sigma reach or exceed
deltaSigma,
                 # linearly interpolate between the nearest points to find
                 # mixed-layer depth
                 if (any(reachedDeltaSigma))
                 {
                     pair <- min(which(reachedDeltaSigma)) + sigmaRefPos - 1
                     linmod <- lm(z[c(pair, pair + 1)] ~ sigma[c(pair,
pair + 1)])
                     mld <- as.vector(linmod$coefficients[1] +
                                      linmod$coefficients[2] *
                                      (sigmaRef + deltaSigma))
                 } else {
                     # otherwise, mixed-layer depth is the deepest point
                     mld <- min(z)
                 }
             } else {
               # if there is no uniform region, just linearly
interpolate mld
               pair <- min(which(diffs >= deltaSigma)) - 1
               linmod <- lm(z[c(pair, pair + 1)] ~ sigma[c(pair, pair + 1)])
               mld <- as.vector(linmod$coefficients[1] +
                                linmod$coefficients[2] *
                                (sigmaRef + deltaSigma))
             }
         } else {
           mld <- min(z)
         }
         if (pos) mld <- abs(mld)
         mld
     }
   }
---------------------------------------
Translate this code into NCL and you are done.

Best fishes for the New Year,
Sam

-- 
email signature Sam McClatchie,
Supervisory oceanographer,
Fisheries oceanography, Ichthyoplankton, Ship Operations
Southwest Fisheries Science Center, NOAA,
8901 La Jolla Shores Dr.
La Jolla, CA 92037-1509, U.S.A.
email <Sam.McClatchie@noaa.gov>
Office: 858 546 7083, Cellular:  858 752 8495
Research home page <www.fishocean.info>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Dec 30 10:31:08 2013