From 6c29188436e3248832492b4f19e7437887a0745f Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Thu, 27 Mar 2025 19:03:52 +1000 Subject: [PATCH] Added geostrophic module --- actions/actioninfo.h | 9 ++- actions/actiontsc.h | 7 +- include/CF.h | 15 ++++ src/CF.cpp | 173 ++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 199 insertions(+), 5 deletions(-) diff --git a/actions/actioninfo.h b/actions/actioninfo.h index b203db5..77a808f 100644 --- a/actions/actioninfo.h +++ b/actions/actioninfo.h @@ -24,10 +24,15 @@ MString ActionInfo_DoAction(const CLArgs& args, D& data) ad = std::move(ret.Value()); } - { // Modules + // Modules + { auto ret = CF::TEOS10Module(args, pars, ad); if(!ret) return "Error adding TEOS-10 derived variables"; } + { + auto ret = CF::GeostrophicModule(args, pars, ad); + if(!ret) return "Error adding geostrophic velocity"; + } const MString titCol = F(C::BWHITE); const MString timeCol = F(C::YELLOW); @@ -72,7 +77,7 @@ MString ActionInfo_DoAction(const CLArgs& args, D& data) // Metainfo message(ad.Is2D(name) ? v2dCol : v3dCol, name, C(), info->lname.Exist() ? (" - " + info->lname) : "", info->targetunit.Exist() ? (", measured in " + unitCol + info->targetunit + C()) : ""); - if(info->comment.Exist()) message(comCol," ", info->comment, C()); + if(info->comment.Exist()) message(comCol, " ", info->comment, C()); if(info->stname.Exist()) message(" Standard name: " + info->stname); // Vertical grid description diff --git a/actions/actiontsc.h b/actions/actiontsc.h index 2f5892d..65f4a48 100644 --- a/actions/actiontsc.h +++ b/actions/actiontsc.h @@ -29,10 +29,15 @@ MString ActionTSC_DoAction(const CLArgs& args, D& ds) ad = std::move(ret.Value()); } - { // Modules + // Modules + { auto ret = CF::TEOS10Module(args, pars, ad); if(!ret) return "Error adding TEOS-10 derived variables"; } + { + auto ret = CF::GeostrophicModule(args, pars, ad); + if(!ret) return "Error adding geostrophic velocity"; + } MString name = args.contains("out") ? args.at("out") : "out.nc"; int compress = args.contains("compress") ? args.at("compress").ToInt() : 3; diff --git a/include/CF.h b/include/CF.h index 5ab3c5f..cdef005 100644 --- a/include/CF.h +++ b/include/CF.h @@ -28,6 +28,19 @@ class CF virtual ~VInfoTEOS10() override = default; }; + struct VInfoGeostrophic: public Adapter::VInfoBase + { + enum Operation + { + U, + V + }; + MString sshvar; + Operation op; + + virtual ~VInfoGeostrophic() override = default; + }; + // Get time variable name static MString GetTVarName(const NCZarr& nc); @@ -94,6 +107,8 @@ class CF } static RetVal> TEOS103DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr proj, std::shared_ptr vert, size_t it); + static RetVal> Geostrophic2DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr proj, size_t it); static Error TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad); + static Error GeostrophicModule(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad); }; diff --git a/src/CF.cpp b/src/CF.cpp index 2792b78..7d0192f 100644 --- a/src/CF.cpp +++ b/src/CF.cpp @@ -456,7 +456,7 @@ Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListE RetVal> CF::TEOS103DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr proj, std::shared_ptr vert, size_t it) { - static const MString pref = "CF::TempSalDer3DReader"; + static const MString pref = "CF::Teos103DReader"; const auto* v = dynamic_cast(vinfo); if(v == nullptr) return {pref, "Incorrect argument type"}; @@ -550,9 +550,119 @@ RetVal> CF::TEOS103DReader(Adapter& ad, const Adapter::V return out; } +RetVal> CF::Geostrophic2DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr proj, size_t it) +{ + static const MString pref = "CF::Geostrophic2DReader"; + + const auto* v = dynamic_cast(vinfo); + if(v == nullptr) return {pref, "Incorrect argument type"}; + + std::shared_ptr out(new Data2D(proj, {.unit = v->targetunit, .stname = v->stname, .lname = v->lname, .comment = v->comment})); + + std::shared_ptr ssh; + { + auto retssh = ad.Read2D(v->sshvar, ad.CurrentTime()); + if(!retssh) return {pref, "Can't read ssh variable " + v->sshvar}; + ssh = retssh.Value(); + } + + if(ssh->Ny() != out->Ny()) return {pref, MString("For ssh variable Ny=") + ssh->Ny() + ", for output variable Ny=" + out->Ny()}; + if(ssh->Nx() != out->Nx()) return {pref, MString("For ssh variable Nx=") + ssh->Nx() + ", for output variable Nx=" + out->Nx()}; + + const UtUnit infr(ssh->Unit()), into("cm"); + const UtUnit outfr("cm/s"), outto(out->Unit()); + const auto cnvin = Converter(infr, into); + const auto cnvout = Converter(outfr, outto); + + if(!cnvin) return {pref, "Can't convert units from \"" + ssh->Unit() + "\" to \"cm\""}; + if(!cnvout) return {pref, "Can't convert units from \"cm/s\" to \"" + ssh->Unit() + "\""}; + + const bool needsconvertin = infr != into; + const bool needsconvertout = outfr != outto; + + auto In = [needsconvertin, &cnvin = cnvin](real r) -> real { return needsconvertin ? cnvin(r) : r; }; + auto Out = [needsconvertout, &cnvout = cnvout](real r) -> real { return needsconvertout ? cnvout(r) : r; }; + + switch(v->op) + { + case(VInfoGeostrophic::U): + { + for(size_t iy = 0; iy < ssh->Ny(); iy++) + { + const real phi = michlib::M_PI * ssh->Lat(iy) / 180.0; + const real sphi = michlib::Sin(phi); + const real s2phi = michlib::Sin(2.0 * phi); + const real g = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi); // g in cm/s + const real f = 2.0 * 7.29211 * sphi; // Coriolis parameter multiplied by 100'000 + const real db = (iy == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh->Lat(iy) - ssh->Lat(iy - 1)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter + const real du = + (iy == ssh->Ny() - 1) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh->Lat(iy + 1) - ssh->Lat(iy)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter + + for(size_t ix = 0; ix < ssh->Nx(); ix++) + { + const bool bok = iy != 0 && !ssh->IsFill(ix, iy - 1) && !ssh->IsFill(ix, iy); + const bool uok = iy != ssh->Ny() - 1 && !ssh->IsFill(ix, iy + 1) && !ssh->IsFill(ix, iy); + if(!(bok || uok)) + (*out)(ix, iy) = out->Fillval(); + else if(!bok) + (*out)(ix, iy) = Out(-g * (In((*ssh)(ix, iy + 1)) - In((*ssh)(ix, iy))) / (f * du)); + else if(!uok) + (*out)(ix, iy) = Out(-g * (In((*ssh)(ix, iy)) - In((*ssh)(ix, iy - 1))) / (f * db)); + else + { + const real bder = -g * (In((*ssh)(ix, iy)) - In((*ssh)(ix, iy - 1))) / (f * db); + const real uder = -g * (In((*ssh)(ix, iy + 1)) - In((*ssh)(ix, iy))) / (f * du); + (*out)(ix, iy) = Out(0.5 * (bder + uder)); + } + } + } + break; + } + case(VInfoGeostrophic::V): + { + for(size_t iy = 0; iy < ssh->Ny(); iy++) + { + const real phi = michlib::M_PI * ssh->Lat(iy) / 180.0; + const real sphi = michlib::Sin(phi); + const real cphi = michlib::Cos(phi); + const real s2phi = michlib::Sin(2.0 * phi); + const real g = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi); // g in cm/s + const real f = 2.0 * 7.29211 * sphi; // Coriolis parameter multiplied by 100'000 + + for(size_t ix = 0; ix < ssh->Nx(); ix++) + { + const real dl = + (ix == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh->Lon(ix) - ssh->Lon(ix - 1)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter + const real dr = (ix == ssh->Nx() - 1) + ? 0.0 + : (michlib::M_PI * 6371.0 * (ssh->Lon(ix + 1) - ssh->Lon(ix)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter + + const bool lok = ix != 0 && !ssh->IsFill(ix - 1, iy) && !ssh->IsFill(ix, iy); + const bool rok = ix != ssh->Nx() - 1 && !ssh->IsFill(ix + 1, iy) && !ssh->IsFill(ix, iy); + if(!(lok || rok)) + (*out)(ix, iy) = out->Fillval(); + else if(!lok) + (*out)(ix, iy) = Out(g * (In((*ssh)(ix + 1, iy)) - In((*ssh)(ix, iy))) / (f * dr)); + else if(!rok) + (*out)(ix, iy) = Out(g * (In((*ssh)(ix, iy)) - In((*ssh)(ix - 1, iy))) / (f * dl)); + else + { + const real lder = g * (In((*ssh)(ix, iy)) - In((*ssh)(ix - 1, iy))) / (f * dl); + const real rder = g * (In((*ssh)(ix + 1, iy)) - In((*ssh)(ix, iy))) / (f * dr); + (*out)(ix, iy) = Out(0.5 * (lder + rder)); + } + } + } + break; + } + } + + return out; +} + Error CF::TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad) { - static const MString pref = "TempSalDerModule"; + static const MString pref = "TEOS10Module"; MString sal, temp, ptemp; @@ -708,3 +818,62 @@ Error CF::TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pa return {}; } + +Error CF::GeostrophicModule(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad) +{ + static const MString pref = "GeostrophicModule"; + + MString ssh; + + for(const auto& v: ad.Vars()) + { + if(!ad.Is2D(v.first)) continue; + const auto& info = v.second.info; + if(!ssh.Exist() && info->stname == "sea_surface_height_above_geoid") ssh = v.first; + if(!ssh.Exist() && info->stname == "sea_surface_elevation") ssh = v.first; + } + + if(!ssh.Exist()) return {}; + + const auto& svar = ad.Vars().at(ssh); + const auto& sproj = svar.proj; + const auto& rinfo = svar.info->rinfo; + + if(!sproj->IsRectangular()) return {}; // We deal only with rectangular grids for now + + { // U-component + auto vinfo = new VInfoGeostrophic(); + + vinfo->sshvar = ssh; + vinfo->op = VInfoGeostrophic::U; + + vinfo->stname = "surface_geostrophic_eastward_sea_water_velocity"; + vinfo->lname = "Geostrophic eastward velocity"; + vinfo->comment = "Module Geostrophic: Geostrophic eastward velocity calculated from sea surface height (" + ssh + ")"; + vinfo->targetunit = "cm/s"; + vinfo->rinfo = rinfo; + vinfo->fillval = NAN; + + const bool addsuc = ad.Add2DVariable("ugs", Adapter::VInfo(vinfo), sproj, Geostrophic2DReader); + if(!addsuc) return {pref, "Variable with name ugs already exist"}; + } + + { // V-component + auto vinfo = new VInfoGeostrophic(); + + vinfo->sshvar = ssh; + vinfo->op = VInfoGeostrophic::V; + + vinfo->stname = "surface_geostrophic_northward_sea_water_velocity"; + vinfo->lname = "Geostrophic northward velocity"; + vinfo->comment = "Module Geostrophic: Geostrophic northward velocity calculated from sea surface height (" + ssh + ")"; + vinfo->targetunit = "cm/s"; + vinfo->rinfo = rinfo; + vinfo->fillval = NAN; + + const bool addsuc = ad.Add2DVariable("vgs", Adapter::VInfo(vinfo), sproj, Geostrophic2DReader); + if(!addsuc) return {pref, "Variable with name vgs already exist"}; + } + + return {}; +}