diff --git a/.gitignore b/.gitignore index f97f8a8..62e4f99 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,10 @@ SST_NEW5/*.nc SST_1-1-0a/*.tar.bz2 SST_1-1-0a/*.nc +# SST_1-1-1* files +SST_1-1-1/*.tar.bz2 +SST_1-1-1/*.nc + # Executables mkhurrell1 diff --git a/SST_1-1-1/161020.tar.bz2 b/SST_1-1-1/161020.tar.bz2 new file mode 100644 index 0000000..179b39e Binary files /dev/null and b/SST_1-1-1/161020.tar.bz2 differ diff --git a/SST_1-1-1/161020_log.txt b/SST_1-1-1/161020_log.txt new file mode 100644 index 0000000..705de95 --- /dev/null +++ b/SST_1-1-1/161020_log.txt @@ -0,0 +1,27 @@ +2016-10-20 11:32:10 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/ [32387] -> ".listing" [1] +2016-10-20 11:32:11 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201601.gz [198514] -> "oiv2mon.201601.gz" [1] +2016-10-20 11:32:12 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201602.gz [196622] -> "oiv2mon.201602.gz" [1] +2016-10-20 11:32:12 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201603.gz [197203] -> "oiv2mon.201603.gz" [1] +2016-10-20 11:32:13 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201604.gz [197428] -> "oiv2mon.201604.gz" [1] +2016-10-20 11:32:14 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201605.gz [198664] -> "oiv2mon.201605.gz" [1] +2016-10-20 11:32:16 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201606.gz [199118] -> "oiv2mon.201606.gz" [1] +2016-10-20 11:32:16 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201607.gz [200825] -> "oiv2mon.201607.gz" [1] +2016-10-20 11:32:17 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201608.gz [203063] -> "oiv2mon.201608.gz" [1] +2016-10-20 11:32:18 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201609.gz [203213] -> "oiv2mon.201609.gz" [1] +FINISHED --2016-10-20 11:32:18-- +Downloaded: 9 files, 1.7M in 0s (3342796 GB/s) +2016-10-20 11:32:19 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/ [32387] -> ".listing" [1] +2016-10-20 11:32:20 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201501.gz [198070] -> "oiv2mon.201501.gz" [1] +2016-10-20 11:32:21 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201502.gz [195578] -> "oiv2mon.201502.gz" [1] +2016-10-20 11:32:22 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201503.gz [196056] -> "oiv2mon.201503.gz" [1] +2016-10-20 11:32:22 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201504.gz [196924] -> "oiv2mon.201504.gz" [1] +2016-10-20 11:32:23 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201505.gz [197966] -> "oiv2mon.201505.gz" [1] +2016-10-20 11:32:23 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201506.gz [198512] -> "oiv2mon.201506.gz" [1] +2016-10-20 11:32:24 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201507.gz [199384] -> "oiv2mon.201507.gz" [1] +2016-10-20 11:32:25 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201508.gz [200925] -> "oiv2mon.201508.gz" [1] +2016-10-20 11:32:26 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201509.gz [199175] -> "oiv2mon.201509.gz" [1] +2016-10-20 11:32:27 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201510.gz [200042] -> "oiv2mon.201510.gz" [1] +2016-10-20 11:32:28 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201511.gz [199439] -> "oiv2mon.201511.gz" [1] +2016-10-20 11:32:29 URL: ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/oiv2mon.201512.gz [199890] -> "oiv2mon.201512.gz" [1] +FINISHED --2016-10-20 11:32:29-- +Downloaded: 12 files, 2.3M in 0s (4436748 GB/s) diff --git a/SST_1-1-1/SSTICE.Update.unf.csh b/SST_1-1-1/SSTICE.Update.unf.csh new file mode 100755 index 0000000..9b7f8c8 --- /dev/null +++ b/SST_1-1-1/SSTICE.Update.unf.csh @@ -0,0 +1,1074 @@ +#!/bin/csh + +set date=`date +%y%m%d` ; # Set date for dynamic run + +# This creates the files like +# MODEL.OI2.ice.mnly.yyyymm-YYYYMM.unf.nc +# MODEL.OI2.sst.mnly.yyyymm-YYYYMM.unf.nc +# ==============================NCL==================================== + +### setenv NCARG_ROOT /usr/local/bin +setenv NCARG_ROOT /work/durack1/Shared/150219_AMIPForcingData/NCL/ncl_ncarg-6.3.0.Linux_RHEL6.4_x86_64_nodap_gcc447 +setenv PATH /usr/local/bin/:${PATH} ; # Add wrapit77 to PATH + +cat >! main.ncl << "END_MAIN_NCL" +;; Not needed from 6.2 onwarded +;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" +;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" +;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" + +procedure linInterpFlip( X, nPass, NX, NY, endPt) +; source: coast_land_NearNbor_bilin.ncl +; modified to take 0.5 to 359.5 as input +; use loop to minimize memory +local np, y, nMSG, nt, nx, ny, dimX, ntim +begin + nx = NX + ny = NY + + dimX = dimsizes(X) + ntim = dimX(0) + do nt=0,ntim-1 + x = lonFlip( X(nt,:,:) ) + do np=1,nPass + x = linmsg(x, (/endPt, nx/)) ; interpolate in longitude + if (ny.gt.0) then ; interpolate in y + y = x(lon|:,lat|:) ; reorder for interpolation + y = linmsg(y, (/endPt,ny/)) + + x = y(lat|:,lon|:) ; reorder to original + end if + + nx = nx + 10 ; expand the number of pts + if (np.gt.4) then + ny = ny + 3 + end if + end do + X(nt,:,:) = (/ lonFlip(x) /) + end do + + nMSG = num(ismissing(X)) + print("****************************") + print("linInterpFlip: nMSG="+nMSG) + print("****************************") + if (nMSG.ne.0) then + printVarSummary(X) + do nt=0,ntim-1 + nMSG = num(ismissing(X(nt,:,:))) + if (nMSG.ne.0) then + print("nt="+nt+" nMSG="+nMSG ) + end if + end do + exit + end if +end + +external SSTICE "./sstice.so" +external JIMH "./consistent.so" +external NEAR_NEIGHBOR "./coast_land_NearNbor.so" + +begin ; MAIN NCL DRIVER + netCDF = True + DEBUG = False + PRNT_DEBUG = True + PLOT_DEBUG = False + PLOT_TYPE = "ps" + ; main input directory + ;diri = "/ptmp/shea/SSTICE/NEW_SST/" + ;diro = "/ptmp/shea/SSTICE/" + ;dirm = "/cgd/cas/shea/JHURRELL/" ; mask ... not used + ;diri = "/project/cas/shea/SSTICE/SST_NEW/" + ;diro = "/project/cas/shea/SSTICE/" + ;dirm = "/project/cas/shea/SSTICE/" ; mask ... not used + ;diri = "/work/durack1/Shared/150219_AMIPForcingData/SST_NEW/$date/" + + ;###### UPDATE : REQUIRES UPDATING ## + diri = "/work/durack1/Shared/150219_AMIPForcingData/SST_1-1-1/161020/" ; ## UPDATE : REQUIRES UPDATING ## + diro = "/work/durack1/Shared/150219_AMIPForcingData/SST_1-1-1/" ; ## UPDATE : REQUIRES UPDATING ## + dirm = "/work/durack1/Shared/150219_AMIPForcingData/SST_1-1-1/" ; ## UPDATE : REQUIRES UPDATING ## + ;###### + + film = "lstags.onedeg.dat" ; sst mask [not used] + + fils = systemfunc ("cd "+diri+"; ls oiv2mon.* ") ; the ftp's file names + nfils = dimsizes(fils) + +;;filoi = "MODEL.OI2.ice.mnly.200501-200903.unf.nc" +;;filos = "MODEL.OI2.sst.mnly.200501-200903.unf.nc" + sfxc = stringtochar( get_file_suffix(fils(0),0) ) + YYYYMM_START = chartostring(sfxc(1:)) + sfxc = stringtochar( get_file_suffix(fils(nfils-1),0) ) + YYYYMM_LAST = chartostring(sfxc(1:)) + filoi = "MODEL.OI2.ice.mnly."+YYYYMM_START+"-"+YYYYMM_LAST+".unf.nc" + filos = "MODEL.OI2.sst.mnly."+YYYYMM_START+"-"+YYYYMM_LAST+".unf.nc" + print("filoi="+filoi) + print("filos="+filos) + + ntim = nfils + + nfStrt = 0 ; nfStrt = nfils-1 + nfLast = nfils-1 + + if (PLOT_DEBUG) then + wks = gsn_open_wks(PLOT_TYPE,"test") ; open workstation (plot destination) + gsn_define_colormap(wks,"BlGrYeOrReVi200") ; choose colormap + gray = NhlNewColor(wks,0.8,0.8,0.8) ; add gray to colormap + + res = True + res@cnFillOn = True ; turn on color + res@gsnSpreadColors = True ; use full range of colormap + res@gsnSpreadColorStart = 2 ; start at color 2 + res@gsnSpreadColorEnd = -3 ; don't use added gray + res@cnLinesOn = False ; no contour lines + res@cnLineLabelsOn = False ; no contour lines + res@cnInfoLabelOn = False ; turn off cn info label + ;res@cnFillDrawOrder = "PreDraw" ; draw contours before continents + res@gsnMaximize = True ; maximize plot + ;res@mpFillOn = False + res@mpCenterLonF = 200. + + res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels + res@cnMinLevelValF = 0. ; set min contour level + res@cnMaxLevelValF = 32. ; set max contour level + res@cnLevelSpacingF = 2. ; set contour spacing + + RES = res + RES@cnMinLevelValF = 0. ; set min contour level + RES@cnMaxLevelValF = 100. ; set max contour level + RES@cnLevelSpacingF = 5. ; set contour spacing + + nfStrt = nfils-1 + ntim = nfLast-nfStrt+1 + end if + ;print(fils(nfStrt:nfLast)) + + year = new (ntim, "integer") + month = new (ntim, "integer") + + print("ntim="+ntim+" nfStrt="+nfStrt+" nfLast="+nfLast) + + nt = -1 + do nf=nfStrt,nfLast + tmp_c = stringtochar(fils(nf)) + nt = nt + 1 + year(nt) = stringtointeger((/tmp_c( 8:11)/)) + month(nt) = stringtointeger((/tmp_c(12:13)/)) + delete(tmp_c) + end do + +print(year+" "+month) + ntStrt = 0 + ntLast = ntim + + day = new (ntim, "integer") + day = 16 ; default + hour = new (ntim, "integer") + hour = 12 ; default + minute = new (ntim, "integer") + minute = 0 + sec = new (ntim, "double") + sec = 0.0d0 + + datesec= new (ntim, "integer") + delete(datesec@_FillValue) + datesec@units = "current seconds of current date" + datesec!0 = "time" + datesec = 43200 ; default + + date_frac= new (ntim, "double") + delete(date_frac@_FillValue) + date_frac@units = "yyyymmdd.fraction_of_day" + date_frac!0 = "time" + date_frac = 0.5d0 + + i = ind(month.eq.2) ; ignore leap year + if (.not.any(ismissing(i))) then + day(i) = 15 + end if + delete(i) + + i = ind(month.eq.2 .or. month.eq.4 .or. month.eq.6 .or. \ + month.eq.9 .or. month.eq.11 ) + if (.not.any(ismissing(i))) then + hour(i) = 0 + datesec(i) = 0 + date_frac(i) = 0.d0 + end if + delete(i) + + units = "days since 1800-01-01 00:00:00" + time = ut_inv_calendar(year,month,day,hour,minute,sec, units, 0) + time!0 = "time" + delete(time@_FillValue) + time@information = "middle of month" + printVarSummary(time) + + date = ut_calendar(time, -2) + date!0 = "time" + date@units = "yyyymmdd" + + date_frac = (/ date + date_frac /) + + print("nfStrt="+nfStrt) + print("nfLast="+nfLast) + print(fils(nfStrt:nfLast) \ + +" "+year+" "+month+" "+day+" "+hour \ + +" "+date+" "+datesec+" "+date_frac) + + nlat = 180 + mlon = 360 + + lat = latGlobeFo(nlat, "lat", "latitude", "degrees_north") + lon = lonGlobeFo(mlon, "lon", "longitude", "degrees_east") + + smsg = 1.e20 + sice = -1.8 + ; create various arrays to be used + sst = new ( (/ntim,nlat,mlon/) , float, smsg) + sst!0 = "time" + sst!1 = "lat" + sst!2 = "lon" + sst&time = time + sst&lat = lat + sst&lon = lon + ; not really needed + ; done in fortran + lsMask_OI = new ( (/nlat,mlon/) , float) ; 0=land , 1=ocean + lsMask_OI@long_name = "land-sea mask" + delete(lsMask_OI@_FillValue) + lsMask_OI!0 = "lat" + lsMask_OI!1 = "lon" + lsMask_OI&lat = lat + lsMask_OI&lon = lon + + sst@longName = "Sea Surface Temperature" + sst@units = "degC" + sst@info = "sst-ice consistency enforced" + + ice = sst ; copy meta data + ice@longName = "sea-ice concentration" + ice@units = "%" + + SST = sst ; SST will be the data written to grid + ICE = ice ; ICE will be the data written to grid + + SST@ice = sice + +;******************************************************** +; Loop over NCEP-OI files +; For memory reasons do manipulations on a per/file basis +;******************************************************** + + YYYYMM = date/100 + +;******************************************************** +; Empirical relationship +;******************************************************** + ; Emp=> Empirical + iceEmp = fspan(0,0.90,91) + sstEmp = 9.328*(0.729-iceEmp^3) - 1.8 ; sst = f(ice) + sstEmp@long_name = "SST MAX: EMPIRICAL" + sstEmp@units = "C" + sstCrit = 5.0 ; 9.328*0.729 - 1.8 + + print(fils) + + do nf=nfStrt,nfLast ; loop over all files + FNAME = diri+fils(nf) ; 19 NOV 2007 + sfx = get_file_suffix(FNAME,0) + + print("nf="+nf+" sfx="+sfx+" FNAME="+FNAME) + if (sfx.eq.".gz") then + print("SUFFIX .gz encountered: Better to manually gzip -d") + exit + ;system("gzip -d "+FNAME) + end if + + filc = stringtochar(fils(nf) ) + yyyy = stringtointeger((/filc( 8:11)/)) + mm = stringtointeger((/filc(12:13)/)) + mm_s = (/filc(12:13)/) ; string + + yyyymm = yyyy*100 + mm + nt = ind(yyyymm.eq.YYYYMM) + print("nf="+nf+" nt="+nt+" yyyymm="+yyyymm) + + if (ismissing(nt(0))) then + print("nt is missing: nf="+nf+" yyyymm="+yyyymm) + exit + end if +print("DEBUG A:+++++++++++++++++++++++") +print("fils(nf)="+fils(nf)) + ; read NCEP data + ; lsMask_OI added 3/27/2006 + SSTICE::sstice (dirm, diri, film, fils(nf) \ + ,yyyy, mm, nlat, mlon, ice(nt,:,:) \ + ,sst(nt,:,:),lsMask_OI, sst@_FillValue ) +print("DEBUG B:+++++++++++++++++++++++") + +;******************************************************** +; gross error checks [basically, protect from round off error] +;******************************************************** + ice(nt,:,:) = ice(nt,:,:) > 0.0 + ice(nt,:,:) = ice(nt,:,:) < 100.0 ; % + + dim2d = dimsizes(sst(nt,:,:)) + +;******************************************************** +; Jim H's consistency fortran code [minor additions by DJS] +; This code assume ice@units = % +;******************************************************** + JIMH::ssticejh (nlat, mlon, ice(nt,:,:), sst(nt,:,:), sst@_FillValue ) + ;if (DEBUG) then + print("date="+date(nt)+" min(ice)="+min(ice(nt,:,:)) \ + +" max(ice)="+max(ice(nt,:,:)) \ + +" min(sst)="+min(sst(nt,:,:)) \ + +" max(sst)="+max(sst(nt,:,:)) ) + ;end if + +;======================================================== +; For *clarity*, explicitly do each of the +; following data operations separately. +;======================================================== + +;******************************************************** +; Set all sst < -1.8 to -1.8 [JimH did this in fortran] +; Safety check to make sure this happened +;******************************************************** + sst(nt,:,:) = sst(nt,:,:) > sice + +;******************************************************** +; where ice > 90 [%], set corresponding sst to -1.8 +; f90: where(ice > 90.) sst = -1.8 +; Safety check to make sure this happened +;******************************************************** +; sst(nt,:,:) = where (ice(nt,:,:).ge.90, sice, sst(nt,:,:)) +;******************************************************** + + sst1d = ndtooned( sst(nt,:,:) ) + ice1d = ndtooned( ice(nt,:,:) ) + i = ind(ice1d.ge.90) + if (.not.any(ismissing(i))) then + sst1d(i) = sice + end if + sst(nt,:,:) = onedtond(sst1d, dim2d) + delete(i) + delete(sst1d) + delete(ice1d) + +;******************************************************** +; where ice < 15 [%], reset ice to 0.0 +; Safety check to make sure this happened +;******************************************************** +; ice(nt,:,:) = where (ice(nt,:,:).lt.15, 0.0 , ice(nt,:,:)) +;******************************************************** + + sst1d = ndtooned( sst(nt,:,:) ) + ice1d = ndtooned( ice(nt,:,:) ) + i = ind(ice1d.lt.15) + if (.not.any(ismissing(i))) then + ice1d(i) = 0.0 + end if + ice(nt,:,:) = onedtond(ice1d, dim2d) + delete(i) + delete(sst1d) + delete(ice1d) + +;******************************************************** +; Any place where sst>sstCrit set the ice=0 +; Safety check to make sure this happened +;******************************************************** +; ice(nt,:,:) = where (ice(nt,:,:).gt.0 .and. sst(nt,:,:).ge.sstCrit \ +; , 0.0 , ice(nt,:,:)) +;******************************************************** + + ice1d = ndtooned( ice(nt,:,:) ) + sst1d = ndtooned( sst(nt,:,:) ) + i = ind(ice1d.gt.0 .and. sst1d.ge.sstCrit) + if (.not.any(ismissing(i))) then + ice1d(i) = 0.0 + end if + ice(nt,:,:) = onedtond(ice1d, dim2d) + delete(i) + delete(sst1d) + delete(ice1d) + +;******************************************************** +; Under the assumption that SST is more reliable than +; sea-ice concentration, where (15<=ice<=90) and the +; sst exceed the empirically determined max .. adjust the sea-ice. +; Make sure tenths are used. Then put into ice_pc +;******************************************************** + ; empirical formula expects tenths + ice1d = ndtooned( ice(nt,:,:)*0.01 ) ; tenths + sst1d = ndtooned( sst(nt,:,:) ) + + n1590 = num(ice1d.ge.0.15 .and. ice1d.lt.0.90\ + .and. .not.ismissing(sst1d) ) + i = ind(ice1d.ge.0.15 .and. ice1d.lt.0.90\ + .and. .not.ismissing(sst1d) ) + + if (.not.any(ismissing(i))) then + ; alter only following + ni = dimsizes(i) + do n=0,ni-1 + sstmx = 9.328*(0.729-ice1d(i(n))^3) - 1.8 + if (sst1d(i(n)).gt.sstmx) then + ice1d(i(n)) = exp( log(0.729-((sstmx+1.8)/9.328))/3.0 ) ; tenths + end if + end do + + + sstmx1d = sst1d ; exact copy + sstmx1d(i) = 9.328*(0.729-ice1d(i)^3) - 1.8 ; empirical max sst + k = ind(ice1d.ge.0.15 .and. ice1d.lt.0.90 \ + .and. sst1d.gt.sstmx1d) + dimk = dimsizes(k) + if (.not.any(ismissing(k))) then + ; calculate the empirical ice max + ice1d(k) = exp( log(0.729-((sstmx1d(k)+1.8)/9.328))/3.0 ) ; tenths + end if + + delete(k) + delete(sst1d) + delete(sstmx1d) + end if + + ice(nt,:,:) = onedtond(ice1d*100, dim2d) ; return to % + ice(nt,:,:) = ice(nt,:,:) > 0.0 + + delete(i) + delete(ice1d) + +;******************************************************** +; Any place where ice < 15% set to 0.0 +;******************************************************** +; ice(nt,:,:) = where (ice(nt,:,:).lt.15, 0.0 , ice(nt,:,:)) +;******************************************************** + + ice1d = ndtooned( ice(nt,:,:) ) + i = ind(ice1d.lt.15) + if (.not.any(ismissing(i))) then + ice1d(i) = 0.0 ; % + end if + ice(nt,:,:) = onedtond(ice1d, dim2d) + delete(i) + delete(ice1d) + + end do + + N18 = num(sst.lt.sice) + N15 = num(ice.gt.0. .and. ice.lt.15) + N15C = num(ice.gt.0. .and. ice.lt.15 .and. sst.gt.sstCrit) + N1590 = num(ice.ge.15 .and. ice.lt.90 .and. sst.gt.sstEmp(1) ) + N90 = num(ice.ge.90 .and. sst.gt.sice) + print("AFTER: N15="+N15+" N15C="+N15C+" N18="+N18+" N90="+N90+" N1590="+N1590) + + sst@info = "sst-ice consistency enforced" + ice@info = "sst-ice consistency enforced" + + if (PLOT_DEBUG) then + nt = ntim-1 + res@gsnCenterString = "SST Grid: "+date(nt) + plot = gsn_csm_contour_map_ce(wks,sst(nt,:,:), res) + RES@gsnCenterString = "ICE Grid: "+date(nt) + plot = gsn_csm_contour_map_ce(wks,ice(nt,:,:), RES) + + res_OI = True + res_OI@gsnSpreadColors = True + res_OI@gsnSpreadColorEnd = -3 ; don't use added gray + res_OI@cnFillOn = True + res_OI@cnFillMode = "CellFill" + ;res_OI@cnMissingValFillColor= "yellow" + res_OI@mpFillOn = False + res_OI@mpFillDrawOrder = "PostDraw" + res_OI@gsnCenterString = "SST Missing" + res_OI@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels + res_OI@cnMinLevelValF = 0 ; set min contour level + res_OI@cnMaxLevelValF = 2 ; one less than max + res_OI@cnLevelSpacingF = 1 ; set contour spacing + ;res_OI@lbLabelStrings = ispan(1,3,1) + if (isatt(res,"mpCenterLonF")) then + res_OI@mpCenterLonF = res@mpCenterLonF + end if + plot = gsn_csm_contour_map_ce(wks,lsMask_OI, res_OI) + ;print(lat+" "+lsMask_OI(:,{0.5})+" "+lsMask_OI(:,{179.5}) ) + end if + ; Antarctic mask + ; NCEP: 0=land , ocean=1 + latMask1d = ndtooned(conform(lsMask_OI, lat, 0)) + lsMask1d = ndtooned(lsMask_OI) + iMask = ind(latMask1d.le.-60 .and. lsMask1d.eq.0) + ; force certain values + ice(:,{-45:35} ,:) = 0.0 + + ;wcStrt = systemfunc("date") + do nt=0,ntim-1 + x1d = ndtooned(sst(nt,:,:)) + x1d(iMask) = sice ; Antarctic land + sst(nt,:,:) = onedtond(x1d, (/nlat,mlon/) ) + + x1d = ndtooned(ice(nt,:,:)) + x1d(iMask) = 100. ; Antarctic land + ice(nt,:,:) = onedtond(x1d, (/nlat,mlon/) ) + end do + ;wallClockElapseTime(wcStrt, "MASK", 0) + + if (PLOT_DEBUG) then + nt = ntim-1 + res@gsnCenterString = "set -1.8" + plot = gsn_csm_contour_map_ce(wks,sst(nt,:,:), res) + RES@gsnCenterString = "set 100" + plot = gsn_csm_contour_map_ce(wks,ice(nt,:,:), RES) + end if + + ; will be inserted over Siberia + sst_zon= dim_avg_Wrap(sst(:,:,{0:180})) + sst_zon@long_name = "Zonal Mean SST" + ice_zon= dim_avg_Wrap(ice(:,:,{0:180})) + ice_zon@long_name = "Zonal Mean SEAICE" + ; completly bogus [again] + ; set China/Siberia to the zonal average + ; makes code 'converge' faster + if (PRNT_DEBUG) then + nMsgS = num(ismissing(sst)) ; total number of _FillValue + print("SST: start: nMsgS ="+nMsgS) + nMsgI = num(ismissing(ice)) ; total number of _FillValue + print("ICE: start: nMsgI ="+nMsgI) + end if + wcStrt = systemfunc("date") + +;*************************************************** +; Begin Interpolation over land +;*************************************************** + nDx = 2 ; longitude only + sst = linmsg (sst, (/0,nDx/)); linearly interpolate over small distances + ice = linmsg (ice, (/0,nDx/)) + + if (PRNT_DEBUG) then + nMsgS = num(ismissing(sst)) ; total number of _FillValue + print("SST: linmsg: nDX="+nDx+" : nMsgS ="+nMsgS) + nMsgI = num(ismissing(ice)) ; total number of _FillValue + print("ICE: linmsg: nDX="+nDx+" : nMsgI ="+nMsgI) + end if + + if (PLOT_DEBUG) then + nt = ntim-1 + res@gsnCenterString = "linmsg="+nDx + plot = gsn_csm_contour_map_ce(wks,sst(nt,:,:), res) + plot = gsn_csm_contour_map_ce(wks,ice(nt,:,:), RES) + end if + ; assumes -180 to +180 + nPtx = 3 ; small inland nearest neighbor + nPty = 1 + ; coastnn works better when 179.5W -> 179.5E + ;NEAR_NEIGHBOR::coastnn(sst,mlon,nlat,ntim,sst@_FillValue,nPtx,nPty) + ;NEAR_NEIGHBOR::coastnn(ice,mlon,nlat,ntim,ice@_FillValue,nPtx,nPty) + + do nt=0,ntim-1 ; use loop to minimize memory + tmp = lonFlip(sst(nt:nt,:,:)) ; make NCEP 179.5W -> 179.5E + NEAR_NEIGHBOR::coastnn(tmp,mlon,nlat, 1 ,tmp@_FillValue,nPtx,nPty) + sst(nt:nt,:,:) = (/ lonFlip(tmp) /) + + tmp = lonFlip(ice(nt:nt,:,:)) + NEAR_NEIGHBOR::coastnn(tmp,mlon,nlat, 1 ,tmp@_FillValue,nPtx,nPty) + ice(nt:nt,:,:) = (/ lonFlip(tmp) /) + end do + + if (PRNT_DEBUG) then + nMsgS = num(ismissing(sst)) ; total number of _FillValue + print("SST: NEAR NEIGHBOR: nPtx="+nPtx+" nPty="+nPty+" nMsgS ="+nMsgS) + nMsgI = num(ismissing(ice)) ; total number of _FillValue + print("ICE: NEAR NEIGHBOR: nPtx="+nPtx+" nPty="+nPty+" nMsgI ="+nMsgI) + end if + + if (PLOT_DEBUG) then + nt = ntim-1 + res@gsnCenterString = "NearNeighbor: nPtx="+nPtx+" nPty="+nPty + plot = gsn_csm_contour_map_ce(wks,sst(nt,:,:), res) + plot = gsn_csm_contour_map_ce(wks,ice(nt,:,:), RES) + end if + ; hasten fill area + sst(:,{25:70},{100}) = (/ sst_zon(:,{25:70}) /) ; arbitrary + sst(:,{50:60},{ 55}) = (/ sst_zon(:,{50:60}) /) ; arbitrary + ice(:,{50:60},{ 55}) = (/ ice_zon(:,{50:60}) /) ; arbitrary + ; large scale iterative bilinear interpolation + nx = 10 ; each is one deg + ny = 3 + endPt = 0 + nPass = 8 + linInterpFlip( sst, nPass, nx, ny, endPt) + if (PRNT_DEBUG) then + nMsgS = num(ismissing(sst)) ; total number of _FillValue + print("SST: nPass="+nPass+" : nMsgS ="+nMsgS) + print(" ") + end if + nPass = 8 + linInterpFlip( ice, nPass, nx, ny, endPt) + + wallClockElapseTime(wcStrt, "NearNeighbor-LinInterp", 0) + + if (PRNT_DEBUG) then + ;nMsgS = num(ismissing(sst)) ; total number of _FillValue + ;print("SST: nPass="+nPass+" : nMsgS ="+nMsgS) + ;print(" ") + nMsgI = num(ismissing(ice)) ; total number of _FillValue + print("ICE: nPass="+nPass+" : nMsgI ="+nMsgI) + print(" ") + end if + + if (PLOT_DEBUG) then + nt = ntim-1 + res@gsnCenterString = "linInterp: nPass="+nPass + plot = gsn_csm_contour_map_ce(wks,sst(nt,:,:), res) + plot = gsn_csm_contour_map_ce(wks,ice(nt,:,:), RES) + end if + + printMinMax(sst , True) + print("nMsg(sst)="+num(ismissing(sst))) + printMinMax(ice , True) + print("nMsg(ice)="+num(ismissing(ice))) + + sst@info = "sst-ice consistency enforced" + ice@info = "sst-ice consistency enforced" + +;******************************************************** +; Force consistency: make sure land interp does not screw up +;******************************************************** + do nt=0,ntim-1 + +;******************************************************** +; Set all sst < -1.8 to -1.8 [JimH did this in fortran] +; Safety check to make sure this happened +;******************************************************** + sst(nt,:,:) = sst(nt,:,:) > sice + +;******************************************************** +; where ice > 90 [%], set corresponding sst to -1.8 +; f90: where(ice > 90.) sst = -1.8 +; Safety check to make sure this happened +;******************************************************** +; sst(nt,:,:) = where (ice(nt,:,:).ge.90, sice, sst(nt,:,:)) +;******************************************************** + + sst1d = ndtooned( sst(nt,:,:) ) + ice1d = ndtooned( ice(nt,:,:) ) + i = ind(ice1d.ge.90) + if (.not.any(ismissing(i))) then + sst1d(i) = sice + end if + sst(nt,:,:) = onedtond(sst1d, dim2d) + delete(i) + delete(sst1d) + delete(ice1d) + +;******************************************************** +; where ice < 15 [%], reset ice to 0.0 +; Safety check to make sure this happened +;******************************************************** +; ice(nt,:,:) = where (ice(nt,:,:).lt.15, 0.0 , ice(nt,:,:)) +;******************************************************** + + sst1d = ndtooned( sst(nt,:,:) ) + ice1d = ndtooned( ice(nt,:,:) ) + i = ind(ice1d.lt.15) + if (.not.any(ismissing(i))) then + ice1d(i) = 0.0 + end if + ice(nt,:,:) = onedtond(ice1d, dim2d) + delete(i) + delete(sst1d) + delete(ice1d) + +;******************************************************** +; Any place where sst>sstCrit set the ice=0 +; Safety check to make sure this happened +;******************************************************** +; ice(nt,:,:) = where (ice(nt,:,:).gt.0 .and. sst(nt,:,:).ge.sstCrit \ +; , 0.0 , ice(nt,:,:)) +;******************************************************** + + ice1d = ndtooned( ice(nt,:,:) ) + sst1d = ndtooned( sst(nt,:,:) ) + i = ind(ice1d.gt.0 .and. sst1d.ge.sstCrit) + if (.not.any(ismissing(i))) then + ice1d(i) = 0.0 + end if + ice(nt,:,:) = onedtond(ice1d, dim2d) + delete(i) + delete(sst1d) + delete(ice1d) + +;******************************************************** +; Under the assumption that SST is more reliable than +; sea-ice concentration, where (15<=ice<=90) and the +; sst exceed the empirically determined max .. adjust the sea-ice. +; Make sure tenths are used. Then put into ice_pc +;******************************************************** + ; empirical formula expects tenths + ice1d = ndtooned( ice(nt,:,:)*0.01 ) + sst1d = ndtooned( sst(nt,:,:) ) + + n1590 = num(ice1d.ge.0.15 .and. ice1d.lt.0.90\ + .and. .not.ismissing(sst1d) ) + i = ind(ice1d.ge.0.15 .and. ice1d.lt.0.90\ + .and. .not.ismissing(sst1d) ) + + if (.not.any(ismissing(i))) then + ; alter only following + ni = dimsizes(i) + do n=0,ni-1 + sstmx = 9.328*(0.729-ice1d(i(n))^3) - 1.8 + if (sst1d(i(n)).gt.sstmx) then + ice1d(i(n)) = exp( log(0.729-((sstmx+1.8)/9.328))/3.0 ) ; tenths + end if + end do + + + sstmx1d = sst1d ; exact copy + sstmx1d(i) = 9.328*(0.729-ice1d(i)^3) - 1.8 ; empirical max sst + k = ind(ice1d.ge.0.15 .and. ice1d.lt.0.90 \ + .and. sst1d.gt.sstmx1d) + dimk = dimsizes(k) + if (.not.any(ismissing(k))) then + ; calculate the empirical ice max + ice1d(k) = exp( log(0.729-((sstmx1d(k)+1.8)/9.328))/3.0 ) ; tenths + end if + + delete(k) + delete(sst1d) + delete(sstmx1d) + end if + + ice(nt,:,:) = onedtond(ice1d*100, dim2d) ; return to % + ice(nt,:,:) = ice(nt,:,:) > 0.0 + + delete(i) + delete(ice1d) + +;******************************************************** +; Any place where ice < 15% set to 0.0 +;******************************************************** +; ice(nt,:,:) = where (ice(nt,:,:).lt.15, 0.0 , ice(nt,:,:)) +;******************************************************** + + ice1d = ndtooned( ice(nt,:,:) ) + i = ind(ice1d.lt.15) + if (.not.any(ismissing(i))) then + ice1d(i) = 0.0 ; % + end if + ice(nt,:,:) = onedtond(ice1d, dim2d) + delete(i) + delete(ice1d) + + end do + + + N18 = num(sst.lt.sice) + N15 = num(ice.gt.0. .and. ice.lt.15) + N15C = num(ice.gt.0. .and. ice.lt.15 .and. sst.gt.sstCrit) + N1590 = num(ice.ge.15 .and. ice.lt.90 .and. sst.gt.sstEmp(1) ) + N90 = num(ice.ge.90 .and. sst.gt.sice) + print("AFTER: N15="+N15+" N15C="+N15C+" N18="+N18+" N90="+N90+" N1590="+N1590) + +;******************************************************** +; Write netCDF +;******************************************************** +if (netCDF) then + + nline = inttochar(10) + + fAtt = True + fAtt@title = "Hurrell Consistent SST with no missing values over land" + fAtt@OI_clim = "1971-2000" + fAtt@NCEP_OI_clim = "1971-2000" + + fAtt@story = nline + \ +" NCL version of JimH sst-ice.consistency.f " + nline + \ +" " + nline + \ +"[a] all sst < -1.8 to -1.8 [sea-ice value] " + nline + \ +"[b] where(ice > 0.9) sst = -1.8 " + nline + \ +"[c] adjust the sst-ice concentration data via Jim Hack formula " + nline + \ +" sstm = 9.328 * (0.729-si**3) - 1.8 [where si is frac sea-ice] " + nline + \ +"[d] Use ad-hoc method to interpolate values over land " + nline + \ +" values near coasts set to nearest neighbor " + nline + \ +" function cssgrid used to interpolat [not nice] " + + fAtt@creator = "Dennis Shea, CGD" + fAtt@creation_date = systemfunc( "date" ) + + ncfile = diro+filos + print (ncfile) + system ("/bin/rm -f " + ncfile) ; remove an pre-file + + ncdf = addfile(ncfile,"c") ; "c"reate the netCDF file + + setfileoption(ncdf,"DefineMode",True) ; EFFICIENCY + + fileattdef( ncdf, fAtt ) + + dimNames = (/ "time", "lon", "lat" /) + dimSizes = (/ -1 , mlon, nlat /) + dimUnlim = (/ True , False, False /) + filedimdef( ncdf, dimNames, dimSizes, dimUnlim ) + + filevardef ( ncdf, "time", typeof(time), getvardims(time) ) + filevarattdef( ncdf, "time", time ) + ; Define 1D variables. + filevardef ( ncdf, "lon", typeof(lon), getvardims(lon) ) + filevarattdef( ncdf, "lon", lon ) + + filevardef ( ncdf, "lat", typeof(lat), getvardims(lat) ) + filevarattdef( ncdf, "lat", lat ) + + filevardef ( ncdf, "date", typeof(date), getvardims(date)) + filevarattdef( ncdf, "date", date ) + + filevardef ( ncdf, "datesec", typeof(datesec), getvardims(datesec) ) + filevarattdef( ncdf, "datesec", datesec ) + + filevardef ( ncdf, "date_frac", typeof(date_frac), getvardims(date_frac) ) + filevarattdef( ncdf, "date_frac", date_frac ) + + filevardef ( ncdf, "SST", typeof(sst), getvardims(sst)) + filevarattdef( ncdf, "SST", sst ) + + setfileoption(ncdf,"DefineMode",False) ; (not really necessary) + + ncdf->time = (/time /) + ncdf->lat = (/lat /) + ncdf->lon = (/lon /) + ncdf->date = (/date /) + ncdf->datesec = (/datesec /) + ncdf->date_frac= (/date_frac /) + ncdf->SST = (/sst /) + + +; ------------- sea-ice netCDF + + gAtt = True + gAtt@title = "Hurrell Consistent ICE with no missing values over land" + + gAtt@story = nline + \ +" NCL version of JimH sst-ice.consistency.f " + nline + \ +" " + nline + \ +"[a] all sst < -1.8 to -1.8 [sea-ice value] " + nline + \ +"[b] where(ice > 0.9) sst = -1.8 " + nline + \ +"[c] adjust the sst-ice concentration data via Jim Hack formula " + nline + \ +" sstm = 9.328 * (0.729-si**3) - 1.8 [where si is frac sea-ice] " + nline + \ +"[d] minor sea-ice adjustments " + + gAtt@creators = "Dennis Shea, CGD" + gAtt@creation_date = systemfunc( "date" ) + + NCFILE = diro+filoi + print (NCFILE) + system ("/bin/rm -f " + NCFILE) ; remove an pre-file + + NCDF = addfile(NCFILE,"c") ; "c"reate the netCDF file + + setfileoption(NCDF,"DefineMode",True) ; EFFICIENCY + + fileattdef( NCDF, gAtt ) + + dimNames = (/ "time", "lon", "lat" /) + dimSizes = (/ -1 , mlon, nlat /) + dimUnlim = (/ True , False, False /) + filedimdef( NCDF, dimNames, dimSizes, dimUnlim ) + + filevardef ( NCDF, "time", typeof(time), getvardims(time) ) + filevarattdef( NCDF, "time", time ) + ; Define 1D variables. + filevardef ( NCDF, "lon", typeof(lon), getvardims(lon) ) + filevarattdef( NCDF, "lon", lon ) + + filevardef ( NCDF, "lat", typeof(lat), getvardims(lat) ) + filevarattdef( NCDF, "lat", lat ) + + filevardef ( NCDF, "date", typeof(date), getvardims(date) ) + filevarattdef( NCDF, "date", date ) + + filevardef ( NCDF, "datesec", typeof(datesec), getvardims(time) ) + filevarattdef( NCDF, "datesec", datesec ) + + filevardef ( NCDF, "date_frac", typeof(date_frac), getvardims(time) ) + filevarattdef( NCDF, "date_frac", date_frac ) + + filevardef ( NCDF, "SEAICE", typeof(ice), getvardims(ice)) + filevarattdef( NCDF, "SEAICE", ice ) + + setfileoption(NCDF,"DefineMode",False) ; (not really necessary) + + NCDF->time = (/time /) + NCDF->lat = (/lat /) + NCDF->lon = (/lon /) + NCDF->date = (/date /) + NCDF->datesec = (/datesec /) + NCDF->date_frac= (/date_frac /) + NCDF->SEAICE = (/ice /) +end if + +end +"END_MAIN_NCL" + +# ===========================FORTRAN=================================== + +cat >! sstice.f << "END_SSTICE" +C NCLFORTSTART + subroutine sstice (dirm, diri, film, fildata + + ,yyyy, mm, nlat, mlon, ice, sst, tagls, zmsg ) + implicit none + character*(*) diri, dirm, film, fildata + integer yyyy, mm, nlat, mlon + real ice(mlon,nlat), sst(mlon,nlat), zmsg + real tagls(mlon,nlat) +C NCLEND +c ************************************************************* +c This section reads in the NCEP OI data [Jim Hurrell read.f] +c ************************************************************* + +c +c This subroutine reads a individual yyyymm reynolds OI SST fields +c +c The geo-location of the SST array elements are: +c SST(1,1) = 0.5E, 89.5S +c SST(1,2) = 0.5E, 88.5S +c SST(2,1) = 1.5E, 89.5S +c SST(360,180) = 359.5E, 89.5N +c +c a land/sea mask should be used to mask out OI SST analyzed values +c not located in the ocean, e.g. data file lstags.onedeg.dat +c land=0 ocean=1 +c +c sst - sea surface temperature array (deg C) +c ice - ice concentration array (%) (0-100, >100 = land or coast) +c iyrst - year of start date of analysis +c imst - month of start date of analysis +c idst - day of start date of analysis +c iyrnd - year of end date of analysis +c imnd - month of end date of analysis +c idnd - day of end date of analysis +c ndays - number of days in analysis (start date thru enddate) +c index - analysis version for reference +c xlon - longitude of center of grid square +c xlat - latitude of center of grid square +c tagls - land/sea tag array (0=land, 1=water) + +c NOTES: +c - land values for sst do not necessarily coincide with land values +c from ice analysis + + integer id, jd, krecl + parameter (id=360,jd=180) + parameter (krecl=id*jd*4) + + character*1 cice(id,jd) + character*6 yyyymm + +c c c real sstnew(id,jd),icenew(id,jd), tagls(id,jd) +c c c real sstnew(id,jd),icenew(id,jd) + real*8 ix +c c c real sstmax(91),icemax(91) + real si, sstm + + integer iyrst,imst,idst,iyrnd,imnd,idnd,ndays,index + integer i,j,ientry, ic + data ientry /0/ + save ientry +c c c save tagls +c upon initial entry: open and read the land / sea mask + +c ************************************************************* +c This section reads in the land-sea mask. +c The OI analysis is done over all ocean areas and +c the Great Lakes. There is no analysis over land. +c The land values are filled by a Cressman interpolation +c to produce a complete grid for possible interpolation to +c other grids. +c ocean: lstag=1 land: lstag=0 +c ************************************************************* +c For the ice fields, the value 122 represents land or coast. +c Note, the ice land mask is a function of the ice analysis +c and may change periodically. +c ************************************************************* + print *, "ENTER SSTICE" + + if (ientry.eq.0) then + print *,"dirm//film=",trim(dirm)//trim(film) + open(50,file=trim(dirm)//trim(film), convert="big_endian", + * form='unformatted',access='direct',recl=krecl) + print *,"Past open statement for mask" + read(50,rec=1) tagls + print *,"Past read statement" + close(50) + ientry = 1 + print *, "=> LAND_SEA MASK READ OK <=" + end if + print *, "AFTER IENTRY:*************************" + +c ************************************************************* +c This section reads in the data [Jim Hurrell read.f] +c ************************************************************* + +c read sst and sea-ice concentration data + + print *,"diri=",diri + print *,"fildata=",fildata + open(10,file=diri//fildata, convert="big_endian", + * form='unformatted') + print *, "AFTER OPEN" + read(10) iyrst,imst,idst,iyrnd,imnd,idnd,ndays,index + print *, "=>fortran: ",iyrst,imst,idst,iyrnd,imnd,idnd,ndays,index + read(10) ((sst(i,j),i=1,id),j=1,jd) + read(10) ((cice(i,j),i=1,id),j=1,jd) + close (10) + + do j = 1, jd + do i = 1, id + ice(i,j) = float(ichar(cice(i,j))) + enddo + enddo + +c set SST and ice to missing over land + + do j = 1, jd + do i = 1, id + if (tagls(i,j).eq.0) sst(i,j) = zmsg + if (ice(i,j).gt.100.) ice(i,j) = zmsg + enddo + enddo + + return + end +"END_SSTICE" + +# consistent.f and coast_land_NearNbor.f same as +# CGD: /fs/cgd/home0/shea/ncld/ncld2/ncld3/hadley + +## Portland Group compiler no longer available +##WRAPIT -d -pg -fPIC sstice.f +##WRAPIT -d -pg -fPIC consistent.f +## WRAPIT -d -pg -fPIC coast_land_NearNbor.f + +WRAPIT -d sstice.f +WRAPIT -d consistent.f +WRAPIT -d coast_land_NearNbor.f + +# ============================Execute=================================== + + ncl main.ncl + +# ============================Clean UP================================== + /bin/rm -f main.ncl # this is local + /bin/rm -f coast_land_NearNbor*o + /bin/rm -f consistent*o + /bin/rm -f objects + /bin/rm -f sstice.f # this is local + /bin/rm -f sstice*o + /bin/rm -f WRAPIT* + /bin/rm -f core + +# ============================copy files================================ + +#scp /ptmp/shea/SSTICE/*nc tramhill.cgd.ucar.edu:/project/cas/shea/hadley/. +exit diff --git a/SST_1-1-1/SST_COMPARE.ncl b/SST_1-1-1/SST_COMPARE.ncl new file mode 100755 index 0000000..288a823 --- /dev/null +++ b/SST_1-1-1/SST_COMPARE.ncl @@ -0,0 +1,93 @@ +; ================================================= +; This script will check the overlay years on +; the update file to make sure the differencess are 0.0 +; For some reason the 1st two months show differences +; over land. No reason why!!!! +; ================================================= +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" + + ; This is the file to be updated + ;dira = "/project/cas/shea/hadley/" + ;dira = "/project/cas/shea/SSTICE/" + dira = "../" + ;fila = "MODEL.SST.HAD187001-198110.OI198111-201103.nc" + ;fila = "MODEL.SST.HAD187001-198110.OI198111-201112.nc" + ;fila = "MODEL.SST.HAD187001-198110.OI198111-201203.nc" + fila = "MODEL.SST.HAD187001-198110.OI198111-201403.nc" + print("fila="+fila) + ; This is the raw update file + ; This includes a few overlap years + ;dirb = "/project/cas/shea/SSTICE/" + dirb = "./" + ;filb = systemfunc("cd "+dirb+" ; ls MODEL.OI2.sst.mnly.*.unf.nc") + ;filb = "MODEL.OI2.sst.mnly.20xxyy-20XXYY.unf.nc" + ;filb = "MODEL.OI2.sst.mnly.201001-201108.unf.nc" + ;filb = "MODEL.OI2.sst.mnly.201001-201203.unf.nc" + ;filb = "MODEL.OI2.sst.mnly.201201-201303.unf.nc" + filb = "MODEL.OI2.sst.mnly.201401-201503.unf.nc" + print("filb="+filb) + + fa = addfile(dira+fila, "r") + fb = addfile(dirb+filb, "r") + + datea = fa->date + dateb = fb->date + + datea = datea/100 ; yyyymm + dateb = dateb/100 + na = dimsizes(datea) + nb = dimsizes(dateb) + + idatea = ind(datea.ge.dateb(0) .and. datea.le.datea(na-1)) + idateb = ind(dateb.ge.dateb(0) .and. dateb.le.datea(na-1)) + print(datea(idatea)+" "+dateb(idateb)) + + ssta = fa->SST(idatea,:,:) + sstb = fb->SST(idateb,:,:) + + printVarSummary(ssta) + printVarSummary(sstb) + + printMinMax(ssta, True) + printMinMax(sstb, True) + print("=================================") + print("=================================") + print("=================================") + + dims = dimsizes(ssta) + print(dims) + ntim = dims(0) + + diff = sstb-ssta + mxdiff = max(abs(diff)) + print("Max diff over ALL dates: ="+mxdiff) + + do nt=0,ntim-1 + print(dateb(idateb(nt)) +" mxdiff="+ max(abs(diff(nt,:,:)))) + end do + + diff@long_name = "SST Diff" + copy_VarCoords(ssta, diff) + +;************************************************ +; create plot +;************************************************ + ;wks = gsn_open_wks("x11","TEST") ; open a ps file + wks = gsn_open_wks("pdf","TEST") ; open a ps file + res = True ; plot mods desired + res@gsnMaximize = True + res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels + res@cnMinLevelValF = -4. ; set min contour level + res@cnMaxLevelValF = 4. ; set max contour level + res@cnLevelSpacingF = 0.5 ; set contour spacing + + ;res@tiMainString = "CCM2 T42 July" ; plot title + ;res@cnInfoLabelOrthogonalPosF = -0.07 ; move the label inside th plot + + do nt=0,ntim-1 + res@gsnCenterString = dateb(idateb(nt)) + plot = gsn_csm_contour_map_ce(wks,diff(nt,:,:), res) ; create plot + end do + diff --git a/SST_1-1-1/coast_land_NearNbor.f b/SST_1-1-1/coast_land_NearNbor.f new file mode 100755 index 0000000..a89be38 --- /dev/null +++ b/SST_1-1-1/coast_land_NearNbor.f @@ -0,0 +1,70 @@ +C NCLFORTSTART + subroutine coastnn (s,mlon,nlat,ntim,smsg,npassx,npassy) + implicit none + integer mlon, nlat, ntim, npassx, npassy + real s(mlon,nlat,ntim), smsg +C NCLEND + real stmp(mlon,nlat) + integer ml, nl, nt, np +c c c integer kmsg, new, nwe + +c physically the idea is that SSTs don't +c . vary much in the east-west direction over small distances +c all this does is set the 1st coastal point to the +c . nearest left/right neighbor + + do nt=1,ntim + do np=1,npassx + + do nl=1,nlat + do ml=1,mlon + stmp(ml,nl) = s(ml,nl,nt) + end do + end do + + do nl=1,nlat +c ! west-to-east [Gulf, Kurishio] + do ml=mlon,2,-1 + if (stmp(ml ,nl).ne.smsg .and.stmp(ml-1,nl).eq.smsg)then + s(ml-1,nl,nt) = s(ml,nl,nt) + end if + end do +c ! east-to-west + do ml=1,mlon-1 + if (stmp(ml,nl) .ne.smsg .and.stmp(ml+1,nl).eq.smsg)then + s(ml+1,nl,nt) = s(ml,nl,nt) + end if + end do + end do + + end do + + do np=1,npassy + + do nl=1,nlat + do ml=1,mlon + stmp(ml,nl) = s(ml,nl,nt) + end do + end do + + do ml=1,mlon +c ! south-to-north + do nl=1,nlat-1 + if (stmp(ml ,nl).ne.smsg .and.stmp(ml,nl+1).eq.smsg)then + s(ml,nl+1,nt) = s(ml,nl,nt) + end if + end do +c ! north-to-south + do nl=nlat,2,-1 + if (stmp(ml ,nl).ne.smsg .and.stmp(ml,nl-1).eq.smsg)then + s(ml,nl-1,nt) = s(ml,nl,nt) + end if + end do + + + end do + end do + end do + + return + end diff --git a/SST_1-1-1/consistent.f b/SST_1-1-1/consistent.f new file mode 100755 index 0000000..c41b796 --- /dev/null +++ b/SST_1-1-1/consistent.f @@ -0,0 +1,132 @@ +C WRAPIT -pg consistent_coast.f + +C NCLFORTSTART + subroutine ssticejh (nlat, mlon, ice, sst, zmsg ) + implicit none + integer nlat, mlon + real ice(mlon,nlat), sst(mlon,nlat), zmsg +C NCLEND + +c NOTES: +c - land values for sst do not necessarily coincide with land values +c from ice analysis + + integer id, jd, krecl + parameter (id=360,jd=180) + parameter (krecl=id*jd*4) + + real sstnew(id,jd),icenew(id,jd) + real*8 ix + real sstmax(91),icemax(91) + real si, sstm + + integer i,j,ic + +c ************************************************************* +c This section modifies sea ice data to be consistent with SST: +c [Jim Hurrell sst-ice.consistency.f] +c ************************************************************* +c +c create a maximum SST allowed for a particular sea-ice concentration +c +c Function created by Jim Hack 2/4/02 +c +c sstmax = 9.328 * (0.729-ice**3) - 1.8 +c +c icemax(1) = 0.00 sstmax(1) = 5.0 +c icemax(16) = 0.15 sstmax(16) = 4.97 ! si cutoff in orginal data +c icemax(90) = 0.89 sstmax(90) = -1.57 +c icemax(91) = 0.90 sstmax(91) = -1.8 + + ic = 0 + do ix = 0.0, 0.90, 0.01 + ic = ic + 1 + icemax(ic) = ix + sstmax(ic) = 9.328 * (0.729-ix**3) - 1.8 + enddo + + do j = 1, jd + do i = 1,id + sstnew(i,j) = sst(i,j) + icenew(i,j) = ice(i,j) + end do + end do + + do j = 1, jd + do i = 1,id + +c (a) first, do not allow sst < -1.8 +c (b) set all values with ice > 90% c to -1.8 +c +c THE FOLLOWING IS DONE IN NCL ... 2006 +c + if (sstnew(i,j).ne.zmsg .and. + + sstnew(i,j).lt.-1.8) sstnew(i,j) = -1.8 + + if (ice(i,j).ne.zmsg.and.sstnew(i,j).ne.zmsg) then + if (ice(i,j).ge.90.0) then + sstnew(i,j) = -1.8 + endif + endif + +c adjust the ice concentration data + + if (sstnew(i,j).ne.zmsg) then + if (icenew(i,j).ne.zmsg.and.icenew(i,j).gt.0.0) then + + if (icenew(i,j).lt.90.0) then ! Don't adjust values > 90% + si = icenew(i,j) * 0.01 ! Convert to fraction from % + sstm = 9.328 * (0.729-si**3) - 1.8 + if (sstnew(i,j).gt.sstm) then + if (sstnew(i,j).gt.sstmax(1)) then + icenew(i,j) = 0.0 + else + do ic = 1, 90 + if (sstnew(i,j).lt.sstmax(ic).and.sstnew(i,j).gt. + * sstmax(ic+1)) then + icenew(i,j) = icemax(ic) * 100. + endif + enddo + endif + endif + endif + + endif + endif + +c force no sea ice < 15%, as in HadISST data +c DJS: UNCOMMENTED 16 Aug 2006 + + if (icenew(i,j).ne.zmsg) then + if (icenew(i,j).lt.15.0) icenew(i,j) = 0.0 +C DJS if (icenew(i,j).gt.99.9) then +C DJS write (*,*) iyr,imn,i,j,ice(i,j) +C DJS endif + endif + enddo + enddo + +c DJS: PUT NEW VALUES BACK TO SST/ICE FOR RETURN TO NCL + + DO J = 1, JD + DO I = 1,ID + SST(I,J) = SSTNEW(I,J) + ICE(I,J) = ICENEW(I,J) + END DO + END DO + +c DJS: THE FOLLOWING IS SOMETHING LIKE WHAT JIM HAD IN hadissst+oiv2.f + +c djs DO J = 1, JD +c djs DO I = 1, ID +c djs IF (ICE(I,J).NE.ZMSG) THEN +c djs IF (ICE(I,J).LT.0.0.OR.ICE(I,J).GT.100.0) THEN +c djs WRITE (*,*) 'OIv2 ',ICE(I,J),I,J +c djs ENDIF +c djs ENDIF +c djs ENDDO +c djs ENDDO + + return + end + diff --git a/SST_1-1-1/download.sh b/SST_1-1-1/download.sh new file mode 100755 index 0000000..c190284 --- /dev/null +++ b/SST_1-1-1/download.sh @@ -0,0 +1,122 @@ +#!/bin/tcsh + +# Author: Paul J. Durack : pauldurack@llnl.gov +# Created on Wed Apr 22 15:42:19 2015 +# @author: durack1 + +# File written to download all OISSTv2 files +# PJD 9 Apr 2015 - Adapted from ERSST_V3b data +# PJD 22 Apr 2015 - Updated to complete data generation +# PJD 23 Apr 2015 - Updated URL and DOI +# PJD 14 Apr 2016 - Updated for V1.0.1 - April 2015 to March 2016 extension +# PJD 25 May 2016 - Updated for V1.0.1 using correct V1.0.0 input data - May 2015 to April 2016 extension +# PJD 26 May 2016 - Added wrapit77 to PATH +# PJD 20 Oct 2016 - Updated for V1.1.1 using V1.1.0 input data - May 2016 to September 2016 extension + +# USER WILL NEED TO SET: +# ENVIRONMENT VARIABLES: NCARG_ROOT (and PATH to include ncl path) +# PATHS: oi2path (and diri, diro and dirm in SSTICE.Update.unf.csh) + +# More info: https://climatedataguide.ucar.edu/climate-data/merged-hadley-noaaoi-sea-surface-temperature-sea-ice-concentration-hurrell-et-al-2008 +# Doc: http://doi.org/10.1175/2008JCLI2292.1 (Hurrell et al., 2008) + +### USER TO SET ### +setenv NCARG_ROOT /work/durack1/Shared/150219_AMIPForcingData/NCL/ncl_ncarg-6.3.0.Linux_RHEL6.4_x86_64_nodap_gcc447 +setenv PATH /work/durack1/Shared/150219_AMIPForcingData/NCL/ncl_ncarg-6.3.0.Linux_RHEL6.4_x86_64_nodap_gcc447/bin:${PATH} +setenv PATH /usr/local/bin/:${PATH} ; # Add wrapit77 to PATH + +###### UPDATE : PATH REQUIRES UPDATING ## +set oi2path=/work/durack1/Shared/150219_AMIPForcingData/SST_1-1-1/ ; ## UPDATE : PATH REQUIRES UPDATING ## +###### + +set date=`date +%y%m%d` +cd ${oi2path} +\rm -r -f ${date} ; # Purge if exists +\mkdir ${date} +cd ${date} + +### Step 1 - get most up-to-date files ### +set currentYear=`date +%y` +set acceptList=oiv2mon.20${currentYear} +echo 'downloading '${acceptList}\*.gz +\wget -o ../${date}_log.txt -nv -nc -nH --cut-dirs=3 -rl1 -A ${acceptList}\*.gz --no-check-certificate ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/ +set previousYear=`expr ${currentYear} - 1` +set previousYear=oiv2mon.20${previousYear} +echo 'downloading '${previousYear}\*.gz +\wget -a ../${date}_log.txt -nv -nc -nH --cut-dirs=3 -rl1 -A ${previousYear}\*.gz --no-check-certificate ftp://ftp.emc.ncep.noaa.gov/cmb/sst/oimonth_v2/ + + +### Step 2 - unzip files ### +\gunzip *.gz + + +### Step 3 - invoke SSTICE.Update.unf.csh ### +# This step requires NCL installed +cd ${oi2path} +###### UPDATE : REQUIRES UPDATING ## +./SSTICE.Update.unf.csh ; ## UPDATE : REQUIRES EDITING TO POINT TO NEW DATA PATH ${date} ## +###### + + +### Step 4 - Interrogate files to make sure things are ok ### +# SST_COMPARE.ncl will generate some slices to peruse + + +### Step 5 - Extract months of data to append to previous files ### +# Extract only 'new' months, Here 12 months and 'time' index values: +# e.g 3,14 (April 2015 through March 2016), 5,15 (June 2015 through April 2016), 16,20 (May 2016 through September 2016) ; ## +###### UPDATE : YEARS REQUIRE UPDATING ## +ncks -O -h -d time,16,20 MODEL.OI2.ice.mnly.201501-201609.unf.nc ICE.update.nc ; # Extract June 2015 to March 2016 ## UPDATE : YEARS IN FILENAMES REQUIRE UPDATING ## +ncks -O -h -d time,16,20 MODEL.OI2.sst.mnly.201501-201609.unf.nc SST.update.nc ; ## UPDATE : YEARS IN FILENAMES REQUIRE UPDATING ## +# Make sure updates went correctly ... only the 12 new months, Check 'date' variable +###### +# Check times of new updates +#ncdump -v date ICE.update.nc +ncdump -v date SST.update.nc +# Check times of older files +#ncdump -v date ../SST_NEW5/MODEL.ICE.HAD187001-198110.OI198111-201604.nc +ncdump -v date ../SST_NEW5/MODEL.SST.HAD187001-198110.OI198111-201604.nc + +### Step 6 - Append new months onto existing data ### +# Purge existing files +#rm -f MODEL.ICE.HAD187001-198110.OI198111-201503.nc ; # Files copied from initial SST_NEW directory - require updating +#rm -f MODEL.SST.HAD187001-198110.OI198111-201503.nc +# Merge old and new files +# V1.0.0 - extend to 201503 +#ncrcat ../MODEL.ICE.HAD187001-198110.OI198111-201403.nc ICE.update.nc MODEL.ICE.HAD187001-198110.OI198111-201503.nc +#ncrcat ../MODEL.SST.HAD187001-198110.OI198111-201403.nc SST.update.nc MODEL.SST.HAD187001-198110.OI198111-201503.nc + + +###### UPDATE : FILENAME/YEARS REQUIRE UPDATING ## +# V1.0.1 - extend to 201604 +ncrcat ../SST_NEW5/MODEL.ICE.HAD187001-198110.OI198111-201604.nc ICE.update.nc MODEL.ICE.HAD187001-198110.OI198111-201609.nc ; ## UPDATE : YEARS IN FILENAMES REQUIRE UPDATING ## +ncrcat ../SST_NEW5/MODEL.SST.HAD187001-198110.OI198111-201604.nc SST.update.nc MODEL.SST.HAD187001-198110.OI198111-201609.nc ; ## UPDATE : YEARS IN FILENAMES REQUIRE UPDATING ## +# Make sure updates went correctly +#ncdump -v date MODEL.ICE.HAD187001-198110.OI198111-201609.nc ; ## UPDATE : YEARS IN FILENAMES REQUIRE UPDATING ## +ncdump -v date MODEL.SST.HAD187001-198110.OI198111-201609.nc ; ## UPDATE : YEARS IN FILENAMES REQUIRE UPDATING ## +# Purge partial new files +rm -f ICE.update.nc ; # These may not purge do to a file handle being unreleased by ncrcat +rm -f MODEL.OI2.ice.mnly.201501-201609.unf.nc ; ## UPDATE : YEARS IN FILENAMES REQUIRE UPDATING ## +rm -f SST.update.nc +rm -f MODEL.OI2.sst.mnly.201501-201609.unf.nc ; ## UPDATE : YEARS IN FILENAMES REQUIRE UPDATING ## +###### + + +### Step 7 - Zip up source data and archive codes ### +echo "Archive source data and purge temp ${date} directory.." +rm -f ${date}.tar.bz2 +tar -cjf ${date}.tar.bz2 ${date} ${date}_log.txt download.sh coast_land_NearNbor.f consistent.f lstags.onedeg.dat SST_COMPARE.ncl SSTICE.Update.unf.csh ; # Archive using bzip2 compression +# Extract using >tar -xjf ${date}.tar.bz2 +# Purge directory +rm -rf ${date} +echo "${date}_AMIP.nc download complete.." +# Conditionally clean up existing versions of files and archive +if ( $1 != "" ) then + \tar -cjf ${1}_archive.tar.bz2 ${1}*.* ; # Archive using bzip2 compression + \rm -f ${1}.tar.bz2 + \rm -f ${1}_log.txt +endif + +: <<-- +comments +-- diff --git a/SST_1-1-1/lstags.onedeg.dat b/SST_1-1-1/lstags.onedeg.dat new file mode 100755 index 0000000..da18367 Binary files /dev/null and b/SST_1-1-1/lstags.onedeg.dat differ