From 5e3c5097c794b6de49f5dbc44d4b785f7aff16c6 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Wed, 6 Mar 2024 13:22:31 -0500 Subject: [PATCH 01/24] support for ERA5 nonaccumulated variables /src/Atmosphere/prepare_era5.jl /src/metadata_era5.jl /examples/forcing_era5.py --- .DS_Store | Bin 0 -> 6148 bytes examples/.DS_Store | Bin 0 -> 6148 bytes examples/forcing_era5.py | 57 +++++++ src/.DS_Store | Bin 0 -> 10244 bytes src/Atmosphere/.DS_Store | Bin 0 -> 6148 bytes src/Atmosphere/prepare_era5.jl | 294 +++++++++++++++++++++++++++++++++ src/metadata_era5.jl | 225 +++++++++++++++++++++++++ 7 files changed, 576 insertions(+) create mode 100644 .DS_Store create mode 100644 examples/.DS_Store create mode 100644 examples/forcing_era5.py create mode 100644 src/.DS_Store create mode 100644 src/Atmosphere/.DS_Store create mode 100644 src/Atmosphere/prepare_era5.jl create mode 100644 src/metadata_era5.jl diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..db0581604fd031ed04616d198bd767ceced4084a GIT binary patch literal 6148 zcmeHKO>fgc5S>j7acGco0I3|1EOD(u(pIVv7dIb=N`O#}-~cGtaRjE0H;SF6D2n7e z{1?vr5cn^g;O$4H*eT*xMb%C;`(}4%JnLtBHcLdTI}bZVO(L?P0O}AcKa8rRIOXLYj(}vvEO^2 z_`=KlY@T-g*;~tMcQszZ zHy9g4oO`n*=E+CQiMSd)qZ8`UfL@SCgCea=WLT%aK`&8iIeiADBg#;8GxRnd8~a^~ ztJ5hm98;g316zm)ili6{7N4gVp6`IyCU>GPeP*vTg*#8IrKK~bk zOwTZ27`RakFstwNySOF2w{F~=*lQi=2dEI?mn)PZX!Le04Y3vPLZx8N5eCrJm@C8x Q#QYJEG?>OPaIFmd1^lw@a{vGU literal 0 HcmV?d00001 diff --git a/examples/.DS_Store b/examples/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..ca00854f28be49e23413595953782a2fad15f730 GIT binary patch literal 6148 zcmeHKO-sZ;41HPmu;Rg91TVt;f(1eEby*Ld{Q>H(A1Kwzf{3@?{W1Q;`jU)`wW|jY zB1;0vn@*E-=1plD05W;H+y@E(BdX$HlV*kJebt@}^30GZmg5!=xWO5&Fz|e&zv9Q5iciIJ)c$e%zPnx#Bpn4TiW>{Uf+cUE^x?9cc{9VaP-Se zy*1(OOQ!7=Yx-D6FIrw1ec(H4Tia8ZU?3O>2G*GY)@+s0mSfmpAQ%V+mJI0kAyE~} zjE$pxI%uo~AnI?V3CG$?XiUmjW^5dJhaw(I^iatyhIlyVQ}!z}HjW++$<2r4Crj>7 z#GcOl$=o5OW7uFI80a#v;g^Nh|Ksob|87uR1p~prS}`EQ>FIRBE#=zU+nmgzW@LL literal 0 HcmV?d00001 diff --git a/examples/forcing_era5.py b/examples/forcing_era5.py new file mode 100644 index 0000000..5b60a3e --- /dev/null +++ b/examples/forcing_era5.py @@ -0,0 +1,57 @@ +# More information about ERA5 products may be found here: +# https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels + +# First, set up the CDS API with these directions: +# https://confluence.ecmwf.int/display/CKB/How+to+install+and+use+CDS+API+on+macOS + +import cdsapi + +c = cdsapi.Client() + +c.retrieve( + 'reanalysis-era5-single-levels', + { + 'product_type': 'reanalysis', + 'format': 'netcdf', + 'variable': [ + '10m_u_component_of_wind', '10m_v_component_of_wind', '2m_dewpoint_temperature', + '2m_temperature', 'mean_eastward_turbulent_surface_stress', 'mean_evaporation_rate', + 'mean_northward_turbulent_surface_stress', 'mean_sea_level_pressure', 'mean_surface_downward_long_wave_radiation_flux', + 'mean_surface_latent_heat_flux', 'mean_surface_net_long_wave_radiation_flux', 'mean_surface_net_short_wave_radiation_flux', + 'mean_surface_sensible_heat_flux', 'mean_total_precipitation_rate', 'total_cloud_cover', + ], + 'year': [ + '2021', + ], + 'month': [ + '08', '09', + ], + 'day': [ + '01', '02', '03', + '04', '05', '06', + '07', '08', '09', + '10', '11', '12', + '13', '14', '15', + '16', '17', '18', + '19', '20', '21', + '22', '23', '24', + '25', '26', '27', + '28', '29', '30', + '31', + ], + 'time': [ + '00:00', '01:00', '02:00', + '03:00', '04:00', '05:00', + '06:00', '07:00', '08:00', + '09:00', '10:00', '11:00', + '12:00', '13:00', '14:00', + '15:00', '16:00', '17:00', + '18:00', '19:00', '20:00', + '21:00', '22:00', '23:00', + ], + 'area': [ + 36, -69, 27, + -62, + ], + }, + 'era5_21.nc') diff --git a/src/.DS_Store b/src/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..90725cc06ffd85b2738c825bbb250dd22affead0 GIT binary patch literal 10244 zcmeI1Pj4GV7{=dq0$JC8WL52f1JX)dI8-HnT2-hQH)$m}Afy^mRV|qHpEz4KyO!-F zZKELH;WO|RIP*c^OCWLQ6rSgujXUd%BMt~fVkX*sckTJ@%skIKlgvg$g59IRmPjNb zAr8h?9VbWO_dG6@O7Y3Zhyi$!sf^`7hBCq7i&wxa;1%!+cm=!y*HHod&SnW~Ci=5i zz$@StSSY~tA;Q7fF5=uoT{>{mEdb&=KGuS3>;sgEir6mV+(bq3OtX6k3N$FT7=gj@ zeZ=aB?IO-iG`Iu?m!N1lD0V1;whk`R>Jqd~^k=VtSKzV&+`He9M>3N`$?g67=V{*W zrFjqEb4oqJaraycu2wL;lj*o%e(jCZWIEg2{5`0wt-tsFMz|5y!f!eU`AKIM&yJ?M z@#F>XzRL5ocX1u}zD|ZG-PWzWJe$QyHk>Gs42G!k@>!A%@{`^CC>uI zJ~%zS_h9Qzv~%y%vpdo0*5(dw9&B%)omIozAAR!hY4`i|Fw1|`?n_XtBhoPUi}O#b z26x<@&~ofn${Fq|G%{7iLS2ex^9b zpj{+=w8@)FW@w4yG`VUro;`8qg1Txyqvd$U63i}bNX zwJuTJ!kPO{{#)>W7$wio;Vk)F^Cm?eFo!P2yyN7N4W;($ApXHDV+NGC8)$Wmo-$`+ zTsb^Z*MG2FZ=&7=+$(ULUuJGPS69sdo2`T7+B(1}Nm6S^o0KJ)<7|hSW7n=@t)cr}^x0~e?M+gt-29n6)uVh&a4+%m~58|hr+@0t}W zxq4@lH_;|Fpwz4qesDKP@lWZ?2Rab01?}*}6JulOc#VnmcYEIFs%A z8*+@E9lKX1b5Lrc+T3@d^IsTk-$@ z>)3VQC$E54;9pW8SbyAkyo-ljyc1h~F5k6BIKIQd!g_NP6~Reg$0O?Nc--)HJl`#{ n?j|D(uG>YNn}|Se*gb}+yDO!D@al2 literal 0 HcmV?d00001 diff --git a/src/Atmosphere/.DS_Store b/src/Atmosphere/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..ba854c0271c5054ca0c82f8dabc953704045a7a9 GIT binary patch literal 6148 zcmeHKyH3ME5S$Guf{-9|G)Vq}NE8%wCJF^L9{?obsRTy`gy>xGKYXH@-91GlHf@Ap zSK6Dq^Uj@a&hmNz$nbS~1oQ#)=!%_fHZ7*l)f?7HGhL$C8rLXsgG>6P5gil1Q31WX zJ{ERw3w*!7Io@!K0Uj|spB0nwtQcD?a7rB1W_F~#7JkXKj@#JsLC%S@>ze)D2K$ct zKDTl4xwi8@UU1L1t;8G!rg+B__nqT`?_cW*a=5kvhhd>WC=d$#I|Z1t)p|RQ(S`z{ zKq#736}T%NFU zbUGvzACh*KRH2CN&h;}3hvbgYh616$ssdX+>}&o%`?>yKCB-Nd2nGI>0@58_jD|c? y%B_vVX|B!aS9CRvs~ndUb}SV$R#Wi--Hh{@c8Gbx%8@g)_z^G}q6q~ys=yZ}C{C;Z literal 0 HcmV?d00001 diff --git a/src/Atmosphere/prepare_era5.jl b/src/Atmosphere/prepare_era5.jl new file mode 100644 index 0000000..6171bd5 --- /dev/null +++ b/src/Atmosphere/prepare_era5.jl @@ -0,0 +1,294 @@ +""" + ROMS.prepare_era5( + atmo_fname,Vnames,filename_prefix,domain_name; + time_origin = DateTime(1858,11,17)) + +Generate ROMS forcing fields from the ERA5 data file `atmo_fname`. + +Translated to Julia by Alexander-Barth and rewritten by lnferris +to use rate variables instead of accumulated variables. + +# Example + +```julia +datadir = "..." +atmo_fname = joinpath(datadir,"ecmwf_sample_data.nc") +filename_prefix = joinpath(datadir,"liguriansea_") +domain_name = "Ligurian Sea Region" +Vnames = ["sustr","svstr","shflux","swflux","swrad","Uwind","Vwind","lwrad", + "lwrad_down","latent","sensible","cloud","rain","Pair","Tair","Qair"] +ROMS.prepare_era5(atmo_fname,Vnames,filename_prefix,domain_name) +) +``` + +Based on forcing/d_ecmwf2roms.m (svn revision 1102): + + Copyright (c) 2002-2017 The ROMS/TOMS Group Hernan G. Arango + Licensed under a MIT/X style license John Wilkin + See License_ROMS.txt + +""" +function prepare_era5( + atmo_fname,Vnames,filename_prefix,domain_name; + time_origin = DateTime(1858,11,17), +) + + +ds_ecmwf = NCDataset(atmo_fname) +lon = ds_ecmwf["longitude"][:] +lat = ds_ecmwf["latitude"][:] +time = ds_ecmwf["time"][:] + +fliplat = lat[2] < lat[1] +if fliplat + lat = reverse(lat) +end + +# The units we have (ECMWFname): https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels +# The units we need (Vname): varinfo.dat + +F = [ + ( + Vname = "sustr", #[N/m2] + output = "sms", + #ECMWFname = "ewss", #[Ns/m2], Eastward turbulent surface stress + #accumulation = true, + #scale = 1.0/Δt_seconds, + ECMWFname = "metss", #[N/m2], Mean eastward turbulent surface stress + accumulation = false, + scale = 1.0, + ), + ( + Vname = "svstr", #[N/m2] + output = "sms", + #ECMWFname = "nsss", #[Ns/m2], Northward turbulent surface stress + #accumulation = true, + #scale = 1.0/Δt_seconds, + ECMWFname = "mntss", #[N/m2], Mean northward turbulent surface stress + accumulation = false, + scale = 1.0, + ), + ( + Vname = "shflux", # [W/m2] + output = "shflux", + #ECMWFname = "", # calculated from Surface sensible heat flux, '' latent heat flux, '' net thermal radiation, '' net solar radiation [J/m2] + #accumulation = true, + #scale = 1.0/Δt_seconds, + ECMWFname = "", # calculated from Mean surface sensible heat flux, '' latent heat flux, '' net long-wave radiation flux, '' net short-wave radiation flux [W/m2] + accumulation = false, + scale = 1.0, + ), + ( + Vname = "swflux", # [m/s] + output = "swflux", + #ECMWFname = "", # calculated from Evaporation [m] and Total precipitation [m] + #accumulation = true, + #scale = 1.0/Δt_seconds, + ECMWFname = "", # calculated from Mean evaporation rate [kg/m2s], Mean total precipitation rate [kg/m2s] + accumulation = false, + scale = 0.001, + ), + ( + Vname = "swrad", # [W/m2] + output = "swrad", + #ECMWFname = "ssr", # [J/m2], Surface net solar radiation + #accumulation = true, + #scale = 1.0/Δt_seconds, + ECMWFname = "msnswrf", # [W/m2], Mean surface net short-wave radiation flux + accumulation = false, + scale = 1.0, + ), + ( + Vname = "Uwind", #[m/s] + output = "wind", + ECMWFname = "u10", #[m/s], 10m u-component of wind + accumulation = false, + scale = 1.0, + ), + ( + Vname = "Vwind", #[m/s] + output = "wind", + ECMWFname = "v10", #[m/s], 10m v-component of wind + accumulation = false, + scale = 1.0, + ), + ( + Vname = "lwrad", # [W/m2] + output = "lwrad", + #ECMWFname = "str", # [J/m2], Surface net thermal radiation + #accumulation = true, + #scale = 1.0/Δt_seconds, + ECMWFname = "msnlwrf", # [W/m2], Mean surface net long-wave radiation flux + accumulation = false, + scale = 1.0, + ), + ( + Vname = "lwrad_down", # [W/m2] + output = "lwrad", + #ECMWFname = "strd", # [J/m2], Surface thermal radiation downwards + #accumulation = true, + #scale = 1.0/Δt_seconds, + ECMWFname = "msdwlwrf", # [W/m2], Mean surface downward long-wave radiation flux + accumulation = false, + scale = 1.0, + ), + ( + Vname = "latent", #[W/m2] + output = "latent", + #ECMWFname = "slhf", #[J/m2], Surface latent heat flux + #accumulation = true, + #scale = 1.0/Δt_seconds, + ECMWFname = "mslhf", #[W/m2], Mean surface latent heat flux + accumulation = false, + scale = 1.0, + ), + ( + Vname = "sensible", #[W/m2] + output = "sensible", + #accumulation = true, + #ECMWFname = "sshf", # [J/m2], Surface sensible heat flux + #scale = 1.0/Δt_seconds, + ECMWFname = "msshf", # [W/m^2], Mean surface sensible heat flux + accumulation = false, + scale = 1.0 + ), + ( + Vname = "cloud", # [fraction] + output = "cloud", + ECMWFname = "tcc", # [fraction], Total cloud cover + accumulation = false, + scale = 1.0, + ), + ( + Vname = "rain", #[kg/m2s] + output = "rain", + #ECMWFname = "tp", #[m], Total precipitation + #scale = 1000.0/Δt_seconds, + #accumulation = true, + ECMWFname = "mtpr", #[kg/m2s], Mean total precipitation rate + accumulation = false, + scale = 1.0, + ), + ( + Vname = "Pair", # [millibar] + output = "Pair", + ECMWFname = "msl", # [Pa], Mean sea level pressure + accumulation = false, + scale = 0.01, + ), + ( + Vname = "Tair", # [degC] + output = "Tair", + ECMWFname = "", # calculated from 2m temperature [K] + accumulation = false, + scale = 1.0, + ), + ( + Vname = "Qair", # [percent] + output = "Qair", + ECMWFname = "", # calculated from 2m temperature [K] and 2m dewpoint temperature [K] + accumulation = false, + scale = 1.0, + ) +] + + +doFields = filter(i -> F[i].Vname in Vnames,1:length(F)) + +filenames = [(F[i].Vname,filename_prefix * "$(F[i].output).nc") for i = doFields] + +# remove existing files +for (Vname,fname) in filenames + if isfile(fname) + rm(fname) + end +end + +for i = doFields + + min_field = Inf + max_field = -Inf + + Vname = F[i].Vname + @info "Processing: $Vname for $(time[1]) - $(time[end])" + Tname = metadata_era5[Vname].Tname + + outfname = filename_prefix * "$(F[i].output).nc" + + ncattrib = OrderedDict( + String(k) => v for (k,v) in + pairs(metadata_era5[Vname].ncattrib)) + + merge!(ncattrib,OrderedDict( + "time" => Tname, + "coordinates" => "lon lat $Tname")) + + ncattrib_time = OrderedDict( + String(k) => v for (k,v) in + pairs(metadata_era5[Tname].ncattrib)) + + merge!(ncattrib_time,OrderedDict( + "units" => "days since $(Dates.format(time_origin,"yyyy-mm-dd HH:MM:SS"))", + "calendar" => "gregorian")) + + dsout = ROMS.def_forcing( + outfname,lon,lat,Vname,Tname,ncattrib,ncattrib_time, + time_origin; + title = "ECMWF Dataset $domain_name from $atmo_fname", + ) + + # Define variables + + dsout["spherical"][:] = 1 # indeed spherical + dsout["lon"][:] = repeat(lon,inner=(1,length(lat))) + dsout["lat"][:] = repeat(lat',inner=(length(lon),1)) + + scale = F[i].scale + previous_field = zeros(length(lon),length(lat)) + + for irec = 1:length(time) + + if Vname == "Tair" + field = nomissing(ds_ecmwf["t2m"][:,:,irec],NaN) + field = field .- 273.15 + elseif Vname == "Qair" + tsur = nomissing(ds_ecmwf["t2m"][:,:,irec],NaN) + tdew = nomissing(ds_ecmwf["d2m"][:,:,irec],NaN) + tsur = tsur .- 273.15 + tdew = tdew .- 273.15 + field = ROMS.relative_humidity.(tsur,tdew) + elseif Vname == "swflux" + evap = Float64.(nomissing(ds_ecmwf["mer"][:,:,irec],NaN)) + prec = Float64.(nomissing(ds_ecmwf["mtpr"][:,:,irec],NaN)) + field = (evap - prec) .* scale; + elseif Vname == "shflux" + sensible = Float64.(nomissing(ds_ecmwf["msshf"][:,:,irec],NaN)) + latent = Float64.(nomissing(ds_ecmwf["mslhf"][:,:,irec],NaN)) + nlwrad = Float64.(nomissing(ds_ecmwf["msnlwrf"][:,:,irec],NaN)) + nsward = Float64.(nomissing(ds_ecmwf["msnswrf"][:,:,irec],NaN)) + field = (sensible + latent + nlwrad + nsward) * F[i].scale + else + field = nomissing(ds_ecmwf[F[i].ECMWFname][:,:,irec],NaN) + field = field * F[i].scale + end + + if fliplat + field = reverse(field,dims=2) + end + + dsout[Tname][irec] = time[irec] + dsout[Vname][:,:,irec] = field + + min_field = min(min_field,minimum(field)) + max_field = max(max_field,maximum(field)) + end + close(dsout) + + @info "Wrote $Vname, Min= $(min_field) Max= $(max_field)" + +end +close(ds_ecmwf) + +return filenames + +end diff --git a/src/metadata_era5.jl b/src/metadata_era5.jl new file mode 100644 index 0000000..f4f384b --- /dev/null +++ b/src/metadata_era5.jl @@ -0,0 +1,225 @@ + +metadata_era5 = Dict( + "cloud_time" => ( + ncattrib = ( + long_name = "cloud fraction time", + ), + ), + "cloud" => ( + Tname = "cloud_time", + ncattrib = ( + long_name = "cloud fraction", + ), + ), + "pair_time" => ( + ncattrib = ( + long_name = "surface air pressure time", + ), + ), + "Pair" => ( + Tname = "pair_time", + ncattrib = ( + long_name = "surface air pressure", + units = "millibar", + ), + ), + "qair_time" => ( + ncattrib = ( + #long_name = "surface air humidity time", + long_name = "2m air humidity time", + ), + ), + "Qair" => ( + Tname = "qair_time", + ncattrib = ( + #long_name = "surface air humidity", + long_name = "2m air humidity", + units = "percentage", + ), + ), + "tair_time" => ( + ncattrib = ( + #long_name = "surface air temperature time", + long_name = "2m air temperature time", + ), + ), + "Tair" => ( + Tname = "tair_time", + ncattrib = ( + #long_name = "surface air temperature", + long_name = "2m air temperature", + units = "Celsius", + ), + ), + "lhf_time" => ( + ncattrib = ( + long_name = "latent heat flux time", + ), + ), + "latent" => ( + Tname = "lhf_time", + ncattrib = ( + long_name = "net latent heat flux", + units = "Watt meter-2", + positive_value = "downward flux, heating", + negative_value = "upward flux, cooling", + ), + ), + + #---- + "lrf_time" => ( + ncattrib = ( + long_name = "longwave radiation flux time", + ), + ), + "lwrad" => ( + Tname = "lrf_time", + ncattrib = ( + long_name = "net longwave radiation flux", + units = "Watt meter-2", + positive_value = "downward flux, heating", + negative_value = "upward flux, cooling", + ), + ), + "lwrad_down" => ( + Tname = "lrf_time", + ncattrib = ( + long_name = "downwelling longwave radiation flux", + units = "Watt meter-2", + positive_value = "downward flux, heating", + negative_value = "upward flux, cooling", + ), + ), + "sen_time" => ( + ncattrib = ( + long_name = "sensible heat flux time", + ), + ), + "sensible" => ( + Tname = "sen_time", + ncattrib = ( + long_name = "net sensible heat flux", + units = "Watt meter-2", + positive_value = "downward flux, heating", + negative_value = "upward flux, cooling", + ), + ), + "shf_time" => ( + ncattrib = ( + long_name = "surface net heat flux time", + ), + ), + "shflux" => ( + Tname = "shf_time", + ncattrib = ( + long_name = "surface net heat flux", + units = "Watt meter-2", + positive_value = "downward flux, heating", + negative_value = "upward flux, cooling", + ), + ), + "sms_time" => ( + ncattrib = ( + long_name = "surface momentum stress time", + ), + ), + "sustr" => ( + Tname = "sms_time", + ncattrib = ( + long_name = "surface u-momentum stress", + units = "Newton meter-2", + ), + ), + "svstr" => ( + Tname = "sms_time", + ncattrib = ( + long_name = "surface v-momentum stress", + units = "Newton meter-2", + ), + long_name_time = "surface momentum stress time", + ), + "swf_time" => ( + ncattrib = ( + long_name = "surface net freshwater flux time", + ), + ), + "swflux" => ( + Tname = "swf_time", + ncattrib = ( + long_name = "surface net freshwater flux (E-P)", + units = "m s-1", + positive_value = "net evaporation", + negative_value = "net precipitation", + ), + ), + "srf_time" => ( + ncattrib = ( + long_name = "shortwave radiation flux time", + ), + ), + "swrad" => ( + Tname = "srf_time", + ncattrib = ( + long_name = "solar shortwave radiation flux", + units = "Watt meter-2", + ), + ), + + "rain_time" => ( + ncattrib = ( + long_name = "rain fall time", + ), + ), + "rain" => ( + Tname = "rain_time", + ncattrib = ( + long_name = "rain fall", + units = "kilogram meter-2 second-1", + ), + ), + + "wind_time" => ( + ncattrib = ( + #long_name = "surface wind time", + long_name = "10m wind time", + ), + ), + "Uwind" => ( + Tname = "wind_time", + ncattrib = ( + #long_name = "surface u-wind component", + long_name = "10m u-wind component", + units = "meter second-1", + ), + ), + + "Vwind" => ( + Tname = "wind_time", + ncattrib = ( + #long_name = "surface v-wind component", + long_name = "10m v-wind component", + units = "meter second-1", + ), + ), + + # "sst_time" => ( + # ncattrib = ( + # long_name = "sea surface temperature time", + # ), + # ), + # "SST" => ( + # Tname = "sst_time", + # ncattrib = ( + # long_name = "sea surface temperature", + # units = "Celsius", + # ), + # ), + # "dQdSST" => ( + # Tname = "sst_time", + # ncattrib = ( + # long_name = "surface net heat flux sensitivity to SST", + # units = "Watt meter-2 Celsius-1", + # ), + # ), + +) From 01246c92c4d5b23554dfe0673c5ee3c0d7631c62 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Mon, 11 Mar 2024 16:23:06 -0400 Subject: [PATCH 02/24] Delete .DS_Store --- .DS_Store | Bin 6148 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 .DS_Store diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index db0581604fd031ed04616d198bd767ceced4084a..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKO>fgc5S>j7acGco0I3|1EOD(u(pIVv7dIb=N`O#}-~cGtaRjE0H;SF6D2n7e z{1?vr5cn^g;O$4H*eT*xMb%C;`(}4%JnLtBHcLdTI}bZVO(L?P0O}AcKa8rRIOXLYj(}vvEO^2 z_`=KlY@T-g*;~tMcQszZ zHy9g4oO`n*=E+CQiMSd)qZ8`UfL@SCgCea=WLT%aK`&8iIeiADBg#;8GxRnd8~a^~ ztJ5hm98;g316zm)ili6{7N4gVp6`IyCU>GPeP*vTg*#8IrKK~bk zOwTZ27`RakFstwNySOF2w{F~=*lQi=2dEI?mn)PZX!Le04Y3vPLZx8N5eCrJm@C8x Q#QYJEG?>OPaIFmd1^lw@a{vGU From fef85e2fa82ea762b13a5b5d6f9e5a76db1180f9 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Mon, 11 Mar 2024 16:23:59 -0400 Subject: [PATCH 03/24] Delete examples/.DS_Store --- examples/.DS_Store | Bin 6148 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 examples/.DS_Store diff --git a/examples/.DS_Store b/examples/.DS_Store deleted file mode 100644 index ca00854f28be49e23413595953782a2fad15f730..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKO-sZ;41HPmu;Rg91TVt;f(1eEby*Ld{Q>H(A1Kwzf{3@?{W1Q;`jU)`wW|jY zB1;0vn@*E-=1plD05W;H+y@E(BdX$HlV*kJebt@}^30GZmg5!=xWO5&Fz|e&zv9Q5iciIJ)c$e%zPnx#Bpn4TiW>{Uf+cUE^x?9cc{9VaP-Se zy*1(OOQ!7=Yx-D6FIrw1ec(H4Tia8ZU?3O>2G*GY)@+s0mSfmpAQ%V+mJI0kAyE~} zjE$pxI%uo~AnI?V3CG$?XiUmjW^5dJhaw(I^iatyhIlyVQ}!z}HjW++$<2r4Crj>7 z#GcOl$=o5OW7uFI80a#v;g^Nh|Ksob|87uR1p~prS}`EQ>FIRBE#=zU+nmgzW@LL From e206009a078a5ee6b0632e6fe9af5420848f2749 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Mon, 11 Mar 2024 16:25:03 -0400 Subject: [PATCH 04/24] Delete src/.DS_Store --- src/.DS_Store | Bin 10244 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/.DS_Store diff --git a/src/.DS_Store b/src/.DS_Store deleted file mode 100644 index 90725cc06ffd85b2738c825bbb250dd22affead0..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 10244 zcmeI1Pj4GV7{=dq0$JC8WL52f1JX)dI8-HnT2-hQH)$m}Afy^mRV|qHpEz4KyO!-F zZKELH;WO|RIP*c^OCWLQ6rSgujXUd%BMt~fVkX*sckTJ@%skIKlgvg$g59IRmPjNb zAr8h?9VbWO_dG6@O7Y3Zhyi$!sf^`7hBCq7i&wxa;1%!+cm=!y*HHod&SnW~Ci=5i zz$@StSSY~tA;Q7fF5=uoT{>{mEdb&=KGuS3>;sgEir6mV+(bq3OtX6k3N$FT7=gj@ zeZ=aB?IO-iG`Iu?m!N1lD0V1;whk`R>Jqd~^k=VtSKzV&+`He9M>3N`$?g67=V{*W zrFjqEb4oqJaraycu2wL;lj*o%e(jCZWIEg2{5`0wt-tsFMz|5y!f!eU`AKIM&yJ?M z@#F>XzRL5ocX1u}zD|ZG-PWzWJe$QyHk>Gs42G!k@>!A%@{`^CC>uI zJ~%zS_h9Qzv~%y%vpdo0*5(dw9&B%)omIozAAR!hY4`i|Fw1|`?n_XtBhoPUi}O#b z26x<@&~ofn${Fq|G%{7iLS2ex^9b zpj{+=w8@)FW@w4yG`VUro;`8qg1Txyqvd$U63i}bNX zwJuTJ!kPO{{#)>W7$wio;Vk)F^Cm?eFo!P2yyN7N4W;($ApXHDV+NGC8)$Wmo-$`+ zTsb^Z*MG2FZ=&7=+$(ULUuJGPS69sdo2`T7+B(1}Nm6S^o0KJ)<7|hSW7n=@t)cr}^x0~e?M+gt-29n6)uVh&a4+%m~58|hr+@0t}W zxq4@lH_;|Fpwz4qesDKP@lWZ?2Rab01?}*}6JulOc#VnmcYEIFs%A z8*+@E9lKX1b5Lrc+T3@d^IsTk-$@ z>)3VQC$E54;9pW8SbyAkyo-ljyc1h~F5k6BIKIQd!g_NP6~Reg$0O?Nc--)HJl`#{ n?j|D(uG>YNn}|Se*gb}+yDO!D@al2 From 4fad40f3ff0157ea32777ce7b70ea34d34a491be Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Mon, 11 Mar 2024 16:25:18 -0400 Subject: [PATCH 05/24] Delete src/Atmosphere/.DS_Store --- src/Atmosphere/.DS_Store | Bin 6148 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/Atmosphere/.DS_Store diff --git a/src/Atmosphere/.DS_Store b/src/Atmosphere/.DS_Store deleted file mode 100644 index ba854c0271c5054ca0c82f8dabc953704045a7a9..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKyH3ME5S$Guf{-9|G)Vq}NE8%wCJF^L9{?obsRTy`gy>xGKYXH@-91GlHf@Ap zSK6Dq^Uj@a&hmNz$nbS~1oQ#)=!%_fHZ7*l)f?7HGhL$C8rLXsgG>6P5gil1Q31WX zJ{ERw3w*!7Io@!K0Uj|spB0nwtQcD?a7rB1W_F~#7JkXKj@#JsLC%S@>ze)D2K$ct zKDTl4xwi8@UU1L1t;8G!rg+B__nqT`?_cW*a=5kvhhd>WC=d$#I|Z1t)p|RQ(S`z{ zKq#736}T%NFU zbUGvzACh*KRH2CN&h;}3hvbgYh616$ssdX+>}&o%`?>yKCB-Nd2nGI>0@58_jD|c? y%B_vVX|B!aS9CRvs~ndUb}SV$R#Wi--Hh{@c8Gbx%8@g)_z^G}q6q~ys=yZ}C{C;Z From 350cf84bb07e713e47e13a4877cb5cae09375bb7 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Fri, 7 Jun 2024 23:46:16 -0400 Subject: [PATCH 06/24] Update generate_config.jl Fix Coriolis --- src/generate_config.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/generate_config.jl b/src/generate_config.jl index a91be4e..0c2fc70 100644 --- a/src/generate_config.jl +++ b/src/generate_config.jl @@ -20,10 +20,11 @@ function generate_config(grid_fname,x,y,h,mask,pm,pn,dndx,dmde,opt) angle = zeros(size(h)); - # Earth rotation - omega = 2*pi/(24*60*60); # Coriolis parameter - f = 2*omega * cos.(pi*y/180); + #omega = 2*pi/(24*60*60); + #f = 2*omega * cos.(pi*y/180); + omega = 7.2921150e-5; + f = 2*omega * sind.(y); @debug "create_grid" create_grid(grid_fname,h,f,x,y,mask,angle,pm,pn,dndx,dmde); From fe7d4ca4768777af8f858d9116617387d48c1ff7 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 01:37:08 -0400 Subject: [PATCH 07/24] x,y scale --- src/create_grid.jl | 27 +++++++-------- src/create_grid_old.jl | 78 ++++++++++++++++++++++++++++++++++++++++++ src/projections.jl | 23 +++++++++++++ src/projections_old.jl | 12 +++++++ 4 files changed, 126 insertions(+), 14 deletions(-) create mode 100644 src/create_grid_old.jl create mode 100644 src/projections_old.jl diff --git a/src/create_grid.jl b/src/create_grid.jl index 7e286d3..d011d68 100644 --- a/src/create_grid.jl +++ b/src/create_grid.jl @@ -23,12 +23,10 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) lon_u,lon_v,lon_psi = stagger(lon_r) lat_u,lat_v,lat_psi = stagger(lat_r) - m = π/180 * earthradius - - x_rho,y_rho = sg_mercator(lon_r,lat_r) - x_psi,y_psi = sg_mercator(lon_psi,lat_psi) - x_u,y_u = sg_mercator(lon_u,lat_u) - x_v,y_v = sg_mercator(lon_v,lat_v) + x_rho,y_rho = map_to_grid(lon_rho,lat_rho,0,0) + x_u,y_u = map_to_grid(lon_u,lat_u,0.5,0) + x_v,y_v = map_to_grid(lon_v,lat_v,0,0.5) + x_psi,y_psi = map_to_grid(lon_psi,lat_psi,0.5,0.5) ds["spherical"][:] = 'T' ds["lat_rho"][:] = lat_r @@ -61,17 +59,18 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) ds["mask_v"][:] = mask_v ds["mask_psi"][:] = mask_psi - ds["x_rho"][:] = x_rho * m - ds["y_rho"][:] = y_rho * m + ds["x_rho"][:] = x_rho + ds["y_rho"][:] = y_rho + + ds["x_psi"][:] = x_psi + ds["y_psi"][:] = y_psi - ds["x_psi"][:] = x_psi * m - ds["y_psi"][:] = y_psi * m + ds["x_u"][:] = x_u + ds["y_u"][:] = y_u - ds["x_u"][:] = x_u * m - ds["y_u"][:] = y_u * m + ds["x_v"][:] = x_v + ds["y_v"][:] = y_v - ds["x_v"][:] = x_v * m - ds["y_v"][:] = y_v * m @debug "closing" close(ds) diff --git a/src/create_grid_old.jl b/src/create_grid_old.jl new file mode 100644 index 0000000..7e286d3 --- /dev/null +++ b/src/create_grid_old.jl @@ -0,0 +1,78 @@ +""" + ROMS.create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) + +Create a NetCDF grid file `fname` using the bathymetry `h`, Coriolis parameter +`f` and longitude, latitude, mask, angle and strechting factors are rho-points. + +!!! note + This function currently only works for non-rotated grids (angle = 0) and + the spherical grids. +""" +function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) + @assert all(angle .== 0) + + xi_rho,eta_rho = size(h) + + if isfile(fname) + rm(fname) + end + + ds = def_grid(fname,xi_rho,eta_rho) + + mask_u,mask_v,mask_psi = stagger_mask(mask_r) + lon_u,lon_v,lon_psi = stagger(lon_r) + lat_u,lat_v,lat_psi = stagger(lat_r) + + m = π/180 * earthradius + + x_rho,y_rho = sg_mercator(lon_r,lat_r) + x_psi,y_psi = sg_mercator(lon_psi,lat_psi) + x_u,y_u = sg_mercator(lon_u,lat_u) + x_v,y_v = sg_mercator(lon_v,lat_v) + + ds["spherical"][:] = 'T' + ds["lat_rho"][:] = lat_r + ds["lat_u"][:] = lat_u + ds["lat_v"][:] = lat_v + ds["lat_psi"][:] = lat_psi + + ds["lon_rho"][:] = lon_r + ds["lon_u"][:] = lon_u + ds["lon_v"][:] = lon_v + ds["lon_psi"][:] = lon_psi + + ds["hraw"][:,:,1] = h + ds["h"][:] = h + + ds["angle"][:] = angle + + ds["xl"][:] = sum(1 ./ pm[:,1]) + ds["el"][:] = sum(1 ./ pn[1,:]) + + ds["f"][:] = f + ds["pm"][:] = pm + ds["pn"][:] = pn + + ds["dndx"][:] = dndx + ds["dmde"][:] = dmde + + ds["mask_rho"][:] = mask_r + ds["mask_u"][:] = mask_u + ds["mask_v"][:] = mask_v + ds["mask_psi"][:] = mask_psi + + ds["x_rho"][:] = x_rho * m + ds["y_rho"][:] = y_rho * m + + ds["x_psi"][:] = x_psi * m + ds["y_psi"][:] = y_psi * m + + ds["x_u"][:] = x_u * m + ds["y_u"][:] = y_u * m + + ds["x_v"][:] = x_v * m + ds["y_v"][:] = y_v * m + + @debug "closing" + close(ds) +end diff --git a/src/projections.jl b/src/projections.jl index b4fefad..91be492 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -10,3 +10,26 @@ function sg_mercator(lon,lat) end return x,y end + + +function map_to_grid(lon,lat,xshift,yshift) + + dlat_m = 110.574e3 + dlon_m = 111.320e3*cosd(lat) + y = similar(lat) + x = similar(lon) + y = 0 + x = 0 + + for j ∈ 2:size(lon,2) + y[:,j] = dlat_m.*(lat[:,j]-lat[:,j-1])+ y[:,j-1] + end + y = y + yshift.*(y[:,2]-y[:,1]) + + for i ∈ 2:size(lon,1) + x[i,:] = dlon_m[i,:].*(lon[i,:]-lon[i-1,:])+ x[i-1,:]; + end + x = x + xshift.*(x[2,:]-x[1,:]); + + return x,y +end diff --git a/src/projections_old.jl b/src/projections_old.jl new file mode 100644 index 0000000..b4fefad --- /dev/null +++ b/src/projections_old.jl @@ -0,0 +1,12 @@ +function sg_mercator(lon::Number,lat::Number) + return lon*pi/180, log(tand(45 + lat/2)) +end + +function sg_mercator(lon,lat) + x = similar(lon) + y = similar(lon) + for i in eachindex(lon) + x[i],y[i] = sg_mercator(lon[i],lat[i]) + end + return x,y +end From d6deef7c09ec7c1a1088c99857d48f3fba20771f Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 01:40:02 -0400 Subject: [PATCH 08/24] Update create_grid.jl --- src/create_grid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/create_grid.jl b/src/create_grid.jl index d011d68..4117611 100644 --- a/src/create_grid.jl +++ b/src/create_grid.jl @@ -23,7 +23,7 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) lon_u,lon_v,lon_psi = stagger(lon_r) lat_u,lat_v,lat_psi = stagger(lat_r) - x_rho,y_rho = map_to_grid(lon_rho,lat_rho,0,0) + x_rho,y_rho = map_to_grid(lon_r,lat_r,0,0) x_u,y_u = map_to_grid(lon_u,lat_u,0.5,0) x_v,y_v = map_to_grid(lon_v,lat_v,0,0.5) x_psi,y_psi = map_to_grid(lon_psi,lat_psi,0.5,0.5) From 488d64ffe94bc42d0aef6d685168345a0f55055c Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 01:43:11 -0400 Subject: [PATCH 09/24] Update projections.jl --- src/projections.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/projections.jl b/src/projections.jl index 91be492..f0b8014 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -15,7 +15,7 @@ end function map_to_grid(lon,lat,xshift,yshift) dlat_m = 110.574e3 - dlon_m = 111.320e3*cosd(lat) + dlon_m = 111.320e3.*cosd(lat) y = similar(lat) x = similar(lon) y = 0 From 560898a391f3021752dc19ae025ae6bb68d9a391 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 01:46:14 -0400 Subject: [PATCH 10/24] Update projections.jl --- src/projections.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/projections.jl b/src/projections.jl index f0b8014..6fb9d51 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -20,16 +20,16 @@ function map_to_grid(lon,lat,xshift,yshift) x = similar(lon) y = 0 x = 0 - - for j ∈ 2:size(lon,2) - y[:,j] = dlat_m.*(lat[:,j]-lat[:,j-1])+ y[:,j-1] - end - y = y + yshift.*(y[:,2]-y[:,1]) - - for i ∈ 2:size(lon,1) - x[i,:] = dlon_m[i,:].*(lon[i,:]-lon[i-1,:])+ x[i-1,:]; - end - x = x + xshift.*(x[2,:]-x[1,:]); + # + # for j ∈ 2:size(lon,2) + # y[:,j] = dlat_m.*(lat[:,j]-lat[:,j-1])+ y[:,j-1] + # end + # y = y + yshift.*(y[:,2]-y[:,1]) + # + # for i ∈ 2:size(lon,1) + # x[i,:] = dlon_m[i,:].*(lon[i,:]-lon[i-1,:])+ x[i-1,:]; + # end + # x = x + xshift.*(x[2,:]-x[1,:]); return x,y end From f245ed075b0c436eca0978a25798ce760bb4f0dc Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 01:50:45 -0400 Subject: [PATCH 11/24] Update projections.jl --- src/projections.jl | 32 ++++++++++++++------------------ 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/src/projections.jl b/src/projections.jl index 6fb9d51..df9ecb5 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -11,25 +11,21 @@ function sg_mercator(lon,lat) return x,y end +function map_to_grid(lon, lat, xshift, yshift) + dlat_m = 110.574e3 + dlon_m = 111.320e3 * cosd.(lat) -function map_to_grid(lon,lat,xshift,yshift) + y = zeros(size(lon)) + for j = 2:size(lon, 2) + y[:, j] = dlat_m * (lat[:, j] .- lat[:, j-1]) + y[:, j-1] + end + y .+= yshift .* (y[:, 2] .- y[:, 1]) - dlat_m = 110.574e3 - dlon_m = 111.320e3.*cosd(lat) - y = similar(lat) - x = similar(lon) - y = 0 - x = 0 - # - # for j ∈ 2:size(lon,2) - # y[:,j] = dlat_m.*(lat[:,j]-lat[:,j-1])+ y[:,j-1] - # end - # y = y + yshift.*(y[:,2]-y[:,1]) - # - # for i ∈ 2:size(lon,1) - # x[i,:] = dlon_m[i,:].*(lon[i,:]-lon[i-1,:])+ x[i-1,:]; - # end - # x = x + xshift.*(x[2,:]-x[1,:]); + x = zeros(size(lon)) + for i = 2:size(lon, 1) + x[i, :] = dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) + x[i-1, :] + end + x .+= xshift .* (x[2, :] .- x[1, :]) - return x,y + return x, y end From 52061046cddefa51cf4cd96ee305e953100d8191 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 01:53:18 -0400 Subject: [PATCH 12/24] Update projections.jl --- src/projections.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/projections.jl b/src/projections.jl index df9ecb5..0073d65 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -17,13 +17,13 @@ function map_to_grid(lon, lat, xshift, yshift) y = zeros(size(lon)) for j = 2:size(lon, 2) - y[:, j] = dlat_m * (lat[:, j] .- lat[:, j-1]) + y[:, j-1] + y[:, j] = dlat_m * (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] end y .+= yshift .* (y[:, 2] .- y[:, 1]) x = zeros(size(lon)) for i = 2:size(lon, 1) - x[i, :] = dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) + x[i-1, :] + x[i, :] = dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) .+ x[i-1, :] end x .+= xshift .* (x[2, :] .- x[1, :]) From 34a45fcd4a8157e4712e7c276d60f7e887e9079f Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 01:56:40 -0400 Subject: [PATCH 13/24] Update projections.jl --- src/projections.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/projections.jl b/src/projections.jl index 0073d65..c5ce2f3 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -15,13 +15,13 @@ function map_to_grid(lon, lat, xshift, yshift) dlat_m = 110.574e3 dlon_m = 111.320e3 * cosd.(lat) - y = zeros(size(lon)) + y = zeros(size(lat)) # Initialize with the same size as lat for j = 2:size(lon, 2) y[:, j] = dlat_m * (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] end y .+= yshift .* (y[:, 2] .- y[:, 1]) - x = zeros(size(lon)) + x = zeros(size(lon)) # Initialize with the same size as lon for i = 2:size(lon, 1) x[i, :] = dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) .+ x[i-1, :] end From 1d9c0b4e19bd6e91acdfd5ff7e1e3eec1208325f Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 02:06:14 -0400 Subject: [PATCH 14/24] Update projections.jl --- src/projections.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/projections.jl b/src/projections.jl index c5ce2f3..3e38e28 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -13,15 +13,21 @@ end function map_to_grid(lon, lat, xshift, yshift) dlat_m = 110.574e3 - dlon_m = 111.320e3 * cosd.(lat) + dlon_m = 111.320e3 .* cosd.(lat) - y = zeros(size(lat)) # Initialize with the same size as lat + # Initialize y with zeros and preallocate memory + y = zeros(size(lat)) + + # Calculate y values for j = 2:size(lon, 2) - y[:, j] = dlat_m * (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] + y[:, j] = dlat_m .* (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] end y .+= yshift .* (y[:, 2] .- y[:, 1]) - x = zeros(size(lon)) # Initialize with the same size as lon + # Initialize x with zeros and preallocate memory + x = zeros(size(lon)) + + # Calculate x values for i = 2:size(lon, 1) x[i, :] = dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) .+ x[i-1, :] end From 5654d9b57db66a2d1aca736c069f3fea901595b6 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 02:09:36 -0400 Subject: [PATCH 15/24] Update projections.jl --- src/projections.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/projections.jl b/src/projections.jl index 3e38e28..e1ae3a7 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -16,20 +16,22 @@ function map_to_grid(lon, lat, xshift, yshift) dlon_m = 111.320e3 .* cosd.(lat) # Initialize y with zeros and preallocate memory - y = zeros(size(lat)) + y = similar(lat) + fill!(y, 0) # Calculate y values for j = 2:size(lon, 2) - y[:, j] = dlat_m .* (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] + y[:, j] .= dlat_m .* (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] end y .+= yshift .* (y[:, 2] .- y[:, 1]) # Initialize x with zeros and preallocate memory - x = zeros(size(lon)) + x = similar(lon) + fill!(x, 0) # Calculate x values for i = 2:size(lon, 1) - x[i, :] = dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) .+ x[i-1, :] + x[i, :] .= dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) .+ x[i-1, :] end x .+= xshift .* (x[2, :] .- x[1, :]) From dd64e68683b810b9b4d02b57cc71b07f4eeba2b8 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 02:17:22 -0400 Subject: [PATCH 16/24] Update projections.jl --- src/projections.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/projections.jl b/src/projections.jl index e1ae3a7..289ba97 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -23,7 +23,7 @@ function map_to_grid(lon, lat, xshift, yshift) for j = 2:size(lon, 2) y[:, j] .= dlat_m .* (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] end - y .+= yshift .* (y[:, 2] .- y[:, 1]) + #y .+= yshift .* (y[:, 2] .- y[:, 1]) # Initialize x with zeros and preallocate memory x = similar(lon) @@ -33,7 +33,7 @@ function map_to_grid(lon, lat, xshift, yshift) for i = 2:size(lon, 1) x[i, :] .= dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) .+ x[i-1, :] end - x .+= xshift .* (x[2, :] .- x[1, :]) + #x .+= xshift .* (x[2, :] .- x[1, :]) return x, y end From 50bd907f61ec9b21fc3ffb6bf62bcc9cfa0242b0 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 02:24:06 -0400 Subject: [PATCH 17/24] Update projections.jl --- src/projections.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/projections.jl b/src/projections.jl index 289ba97..194dccf 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -23,7 +23,9 @@ function map_to_grid(lon, lat, xshift, yshift) for j = 2:size(lon, 2) y[:, j] .= dlat_m .* (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] end - #y .+= yshift .* (y[:, 2] .- y[:, 1]) + for i = 2:size(lon, 1) + y[i,:] .+= yshift .* (y[i, 2] .- y[i, 1]) + end # Initialize x with zeros and preallocate memory x = similar(lon) From 314da5b3cbcb2bd3f8c7f835d6781c087e2129f5 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 02:24:51 -0400 Subject: [PATCH 18/24] Update projections.jl --- src/projections.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/projections.jl b/src/projections.jl index 194dccf..3f23481 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -24,7 +24,7 @@ function map_to_grid(lon, lat, xshift, yshift) y[:, j] .= dlat_m .* (lat[:, j] .- lat[:, j-1]) .+ y[:, j-1] end for i = 2:size(lon, 1) - y[i,:] .+= yshift .* (y[i, 2] .- y[i, 1]) + y[i,:] .+= yshift * (y[i, 2] - y[i, 1]) end # Initialize x with zeros and preallocate memory From 9b5048b87316400730406788303aa3ce67552047 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 02:26:51 -0400 Subject: [PATCH 19/24] Update projections.jl --- src/projections.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/projections.jl b/src/projections.jl index 3f23481..da37923 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -35,7 +35,9 @@ function map_to_grid(lon, lat, xshift, yshift) for i = 2:size(lon, 1) x[i, :] .= dlon_m[i, :] .* (lon[i, :] .- lon[i-1, :]) .+ x[i-1, :] end - #x .+= xshift .* (x[2, :] .- x[1, :]) + for j = 2:size(lon, 2) + x[:,j] .+= xshift * (x[2, j] - x[1, j]) + end return x, y end From a274918372bd0a1e52a7b19a797b00582175e3ff Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sat, 8 Jun 2024 21:14:57 -0400 Subject: [PATCH 20/24] x y grid fix --- src/create_grid.jl | 25 +++++++------ src/create_grid_laur.jl | 77 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 13 deletions(-) create mode 100644 src/create_grid_laur.jl diff --git a/src/create_grid.jl b/src/create_grid.jl index 4117611..4f1f5c4 100644 --- a/src/create_grid.jl +++ b/src/create_grid.jl @@ -23,10 +23,10 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) lon_u,lon_v,lon_psi = stagger(lon_r) lat_u,lat_v,lat_psi = stagger(lat_r) - x_rho,y_rho = map_to_grid(lon_r,lat_r,0,0) - x_u,y_u = map_to_grid(lon_u,lat_u,0.5,0) - x_v,y_v = map_to_grid(lon_v,lat_v,0,0.5) - x_psi,y_psi = map_to_grid(lon_psi,lat_psi,0.5,0.5) + x_rho,y_rho = sg_mercator(lon_r,lat_r) + x_psi,y_psi = sg_mercator(lon_psi,lat_psi) + x_u,y_u = sg_mercator(lon_u,lat_u) + x_v,y_v = sg_mercator(lon_v,lat_v) ds["spherical"][:] = 'T' ds["lat_rho"][:] = lat_r @@ -59,18 +59,17 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) ds["mask_v"][:] = mask_v ds["mask_psi"][:] = mask_psi - ds["x_rho"][:] = x_rho - ds["y_rho"][:] = y_rho + ds["x_rho"][:] = x_rho * earthradius + ds["y_rho"][:] = y_rho * earthradius - ds["x_psi"][:] = x_psi - ds["y_psi"][:] = y_psi + ds["x_psi"][:] = x_psi * earthradius + ds["y_psi"][:] = y_psi * earthradius - ds["x_u"][:] = x_u - ds["y_u"][:] = y_u - - ds["x_v"][:] = x_v - ds["y_v"][:] = y_v + ds["x_u"][:] = x_u * earthradius + ds["y_u"][:] = y_u * earthradius + ds["x_v"][:] = x_v * earthradius + ds["y_v"][:] = y_v * earthradius @debug "closing" close(ds) diff --git a/src/create_grid_laur.jl b/src/create_grid_laur.jl new file mode 100644 index 0000000..4117611 --- /dev/null +++ b/src/create_grid_laur.jl @@ -0,0 +1,77 @@ +""" + ROMS.create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) + +Create a NetCDF grid file `fname` using the bathymetry `h`, Coriolis parameter +`f` and longitude, latitude, mask, angle and strechting factors are rho-points. + +!!! note + This function currently only works for non-rotated grids (angle = 0) and + the spherical grids. +""" +function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) + @assert all(angle .== 0) + + xi_rho,eta_rho = size(h) + + if isfile(fname) + rm(fname) + end + + ds = def_grid(fname,xi_rho,eta_rho) + + mask_u,mask_v,mask_psi = stagger_mask(mask_r) + lon_u,lon_v,lon_psi = stagger(lon_r) + lat_u,lat_v,lat_psi = stagger(lat_r) + + x_rho,y_rho = map_to_grid(lon_r,lat_r,0,0) + x_u,y_u = map_to_grid(lon_u,lat_u,0.5,0) + x_v,y_v = map_to_grid(lon_v,lat_v,0,0.5) + x_psi,y_psi = map_to_grid(lon_psi,lat_psi,0.5,0.5) + + ds["spherical"][:] = 'T' + ds["lat_rho"][:] = lat_r + ds["lat_u"][:] = lat_u + ds["lat_v"][:] = lat_v + ds["lat_psi"][:] = lat_psi + + ds["lon_rho"][:] = lon_r + ds["lon_u"][:] = lon_u + ds["lon_v"][:] = lon_v + ds["lon_psi"][:] = lon_psi + + ds["hraw"][:,:,1] = h + ds["h"][:] = h + + ds["angle"][:] = angle + + ds["xl"][:] = sum(1 ./ pm[:,1]) + ds["el"][:] = sum(1 ./ pn[1,:]) + + ds["f"][:] = f + ds["pm"][:] = pm + ds["pn"][:] = pn + + ds["dndx"][:] = dndx + ds["dmde"][:] = dmde + + ds["mask_rho"][:] = mask_r + ds["mask_u"][:] = mask_u + ds["mask_v"][:] = mask_v + ds["mask_psi"][:] = mask_psi + + ds["x_rho"][:] = x_rho + ds["y_rho"][:] = y_rho + + ds["x_psi"][:] = x_psi + ds["y_psi"][:] = y_psi + + ds["x_u"][:] = x_u + ds["y_u"][:] = y_u + + ds["x_v"][:] = x_v + ds["y_v"][:] = y_v + + + @debug "closing" + close(ds) +end From 99d4c38bd83d6ba39ca68f4eed08d2952ac1a668 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sun, 9 Jun 2024 00:07:34 -0400 Subject: [PATCH 21/24] swap --- src/create_grid.jl | 25 ++++++++++--------- ...reate_grid_laur.jl => create_grid_alex.jl} | 25 +++++++++---------- 2 files changed, 25 insertions(+), 25 deletions(-) rename src/{create_grid_laur.jl => create_grid_alex.jl} (74%) diff --git a/src/create_grid.jl b/src/create_grid.jl index 4f1f5c4..4117611 100644 --- a/src/create_grid.jl +++ b/src/create_grid.jl @@ -23,10 +23,10 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) lon_u,lon_v,lon_psi = stagger(lon_r) lat_u,lat_v,lat_psi = stagger(lat_r) - x_rho,y_rho = sg_mercator(lon_r,lat_r) - x_psi,y_psi = sg_mercator(lon_psi,lat_psi) - x_u,y_u = sg_mercator(lon_u,lat_u) - x_v,y_v = sg_mercator(lon_v,lat_v) + x_rho,y_rho = map_to_grid(lon_r,lat_r,0,0) + x_u,y_u = map_to_grid(lon_u,lat_u,0.5,0) + x_v,y_v = map_to_grid(lon_v,lat_v,0,0.5) + x_psi,y_psi = map_to_grid(lon_psi,lat_psi,0.5,0.5) ds["spherical"][:] = 'T' ds["lat_rho"][:] = lat_r @@ -59,17 +59,18 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) ds["mask_v"][:] = mask_v ds["mask_psi"][:] = mask_psi - ds["x_rho"][:] = x_rho * earthradius - ds["y_rho"][:] = y_rho * earthradius + ds["x_rho"][:] = x_rho + ds["y_rho"][:] = y_rho - ds["x_psi"][:] = x_psi * earthradius - ds["y_psi"][:] = y_psi * earthradius + ds["x_psi"][:] = x_psi + ds["y_psi"][:] = y_psi - ds["x_u"][:] = x_u * earthradius - ds["y_u"][:] = y_u * earthradius + ds["x_u"][:] = x_u + ds["y_u"][:] = y_u + + ds["x_v"][:] = x_v + ds["y_v"][:] = y_v - ds["x_v"][:] = x_v * earthradius - ds["y_v"][:] = y_v * earthradius @debug "closing" close(ds) diff --git a/src/create_grid_laur.jl b/src/create_grid_alex.jl similarity index 74% rename from src/create_grid_laur.jl rename to src/create_grid_alex.jl index 4117611..4f1f5c4 100644 --- a/src/create_grid_laur.jl +++ b/src/create_grid_alex.jl @@ -23,10 +23,10 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) lon_u,lon_v,lon_psi = stagger(lon_r) lat_u,lat_v,lat_psi = stagger(lat_r) - x_rho,y_rho = map_to_grid(lon_r,lat_r,0,0) - x_u,y_u = map_to_grid(lon_u,lat_u,0.5,0) - x_v,y_v = map_to_grid(lon_v,lat_v,0,0.5) - x_psi,y_psi = map_to_grid(lon_psi,lat_psi,0.5,0.5) + x_rho,y_rho = sg_mercator(lon_r,lat_r) + x_psi,y_psi = sg_mercator(lon_psi,lat_psi) + x_u,y_u = sg_mercator(lon_u,lat_u) + x_v,y_v = sg_mercator(lon_v,lat_v) ds["spherical"][:] = 'T' ds["lat_rho"][:] = lat_r @@ -59,18 +59,17 @@ function create_grid(fname,h,f,lon_r,lat_r,mask_r,angle,pm,pn,dndx,dmde) ds["mask_v"][:] = mask_v ds["mask_psi"][:] = mask_psi - ds["x_rho"][:] = x_rho - ds["y_rho"][:] = y_rho + ds["x_rho"][:] = x_rho * earthradius + ds["y_rho"][:] = y_rho * earthradius - ds["x_psi"][:] = x_psi - ds["y_psi"][:] = y_psi + ds["x_psi"][:] = x_psi * earthradius + ds["y_psi"][:] = y_psi * earthradius - ds["x_u"][:] = x_u - ds["y_u"][:] = y_u - - ds["x_v"][:] = x_v - ds["y_v"][:] = y_v + ds["x_u"][:] = x_u * earthradius + ds["y_u"][:] = y_u * earthradius + ds["x_v"][:] = x_v * earthradius + ds["y_v"][:] = y_v * earthradius @debug "closing" close(ds) From fe793d1b05e86f822be9798325a4d79ba2b20281 Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sun, 9 Jun 2024 01:13:20 -0400 Subject: [PATCH 22/24] Update ROMS.jl --- src/ROMS.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ROMS.jl b/src/ROMS.jl index c63fe39..148231e 100644 --- a/src/ROMS.jl +++ b/src/ROMS.jl @@ -55,6 +55,7 @@ include("opendap.jl") include("HYCOM.jl") include("Atmosphere/thermodynamics.jl") include("Atmosphere/prepare_ecmwf.jl") +include("Atmosphere/prepare_era5.jl") include("Atmosphere/prepare_gfs.jl") end From 4aa55c82d783cfe8f6a69ce223ff2eb35b8053ab Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sun, 9 Jun 2024 01:15:47 -0400 Subject: [PATCH 23/24] Update ROMS.jl --- src/ROMS.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ROMS.jl b/src/ROMS.jl index 148231e..20704c5 100644 --- a/src/ROMS.jl +++ b/src/ROMS.jl @@ -35,6 +35,7 @@ include("vavg.jl") include("extract_ic.jl") include("read_time.jl") include("metadata.jl") +include("metadata_era5.jl") include("extract_bc.jl") include("generate_grid.jl") include("infile.jl") From cffd8a25826a8dafabe9c7d44cf6b931027fcfdc Mon Sep 17 00:00:00 2001 From: Laur Ferris Date: Sun, 9 Jun 2024 12:54:18 -0400 Subject: [PATCH 24/24] Update projections.jl --- src/projections.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/projections.jl b/src/projections.jl index da37923..f1e5672 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -39,5 +39,10 @@ function map_to_grid(lon, lat, xshift, yshift) x[:,j] .+= xshift * (x[2, j] - x[1, j]) end + # center the x-grid + for j = 1:size(lon,2) + x[:, j] .+= - dlon_m[:, j].* (lon[end, j] - lon[1, j])/2 + end + return x, y end