diff --git a/Aforc_ERA5/era5_crocotools_param.py b/Aforc_ERA5/era5_crocotools_param.py index ddc5c9fe7dfff7ba0fbfbd6e277df34b6803ea03..fd3630ea0b6b8ec8467fba1565aee1f29ba3d4aa 100644 --- a/Aforc_ERA5/era5_crocotools_param.py +++ b/Aforc_ERA5/era5_crocotools_param.py @@ -8,7 +8,6 @@ # # General path # -##config_dir = '../croco/Run_TEST/' # must be the same than crocotools_param config_dir = 'path_to_my_run_dir/' # must be the same than crocotools_param config_name = 'Benguela_LR' # diff --git a/Aforc_ERA5/make_bry_wkb_ERA5.m b/Aforc_ERA5/make_bry_wkb_ERA5.m new file mode 100644 index 0000000000000000000000000000000000000000..b9c163c899e4dd9d93833f5b656b65e5be4b51d6 --- /dev/null +++ b/Aforc_ERA5/make_bry_wkb_ERA5.m @@ -0,0 +1,190 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Build a CROCO boundary file +% +% Extrapole and interpole temperature and salinity from a +% climatology to get boundary conditions for +% CROCO (boundary netcdf file) . +% Get the velocities and sea surface elevation via a +% geostrophic computation. +% +% Data input format (netcdf): +% temperature(T, Z, Y, X) +% T : time [Months] +% Z : Depth [m] +% Y : Latitude [degree north] +% X : Longitude [degree east] +% +% Data source : IRI/LDEO climate Data Library (World Ocean Atlas 1998) +% http://ingrid.ldgo.columbia.edu/ +% http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NODC/.WOA98/ +% +% Further Information: +% http://www.croco-ocean.org +% +% This file is part of CROCOTOOLS +% +% CROCOTOOLS is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published +% by the Free Software Foundation; either version 2 of the License, +% or (at your option) any later version. +% +% CROCOTOOLS is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place, Suite 330, Boston, +% MA 02111-1307 USA +% +% Copyright (c) 2005-2006 by Pierrick Penven +% e-mail:Pierrick.Penven@ird.fr +% +% Updated 2016 by Patrick Marchesiello & Rachid Benshila +% for wave forcing +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +clear all +close all +%%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%% +% +% Common parameters +% +crocotools_param +% +%%%%%%%%%%%%%%%%%%% END USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%% +disp(' ') +disp([' Title: ',CROCO_title]) +% +% Read in the grid +% +disp(' ') +disp(' Read in the grid...') +nc=netcdf(grdname); +h=nc{'h'}(:); +lon=nc{'lon_rho'}(:); +lat=nc{'lat_rho'}(:); +angle=nc{'angle'}(:); +mask=nc{'mask_rho'}(:); +Lp=length(nc('xi_rho')); +Mp=length(nc('eta_rho')); +hmax=max(max(nc{'h'}(:))); +result=close(nc); +grid_angle=mean(mean(angle(mask==1))); +[M L]=size(h); +% +% Loop over monthly files +% +for Y=Ymin:Ymax + if Y==Ymin + mo_min=Mmin; + else + mo_min=1; + end + if Y==Ymax + mo_max=Mmax; + else + mo_max=12; + end + for M=mo_min:mo_max +% +% Forcing file name +% + frc_prefix=[frc_prefix,'_ERA5_']; + if level==0 + nc_suffix='.nc'; + else + nc_suffix=['.nc.',num2str(level)]; + end + frcname=[frc_prefix,'Y',num2str(Y),... + 'M',num2str(M),nc_suffix]; +% +% WKB file name +% + wkb_prefix=[wkb_prefix,'_ERA5_']; + brywkbname=[wkb_prefix,'Y',num2str(Y),... + 'M',num2str(M),nc_suffix]; + disp(' ') + disp([' Making file: ',brywkbname]) + disp([' from: ',frcname]) +% +% Read wave fields data +% + nc=netcdf(frcname); + Awave=nc{'Awave'}(:); + Dwave=nc{'Dwave'}(:); + Pwave=nc{'Pwave'}(:); + time =nc{'wwv_time'}(:); + close(nc) + + wkb_time=time; + wkb_cycle=wkb_time(end); +% +% Create the boundary file +% + create_bryfile_wkb(brywkbname,grdname,CROCO_title,wkb_obc,... + wkb_time,wkb_cycle,'clobber'); + disp(' ') +% +% +% Extract boundary data: loop on time and boundaries +% note: in ERA5, meteo convention for wave direction: +% --> zero means "coming from north" and 90 "coming from east" +% + disp(' Extracting and writing file ...') + nc=netcdf(brywkbname,'write'); + g=9.81; + for l=1:length(time) + for obcndx=1:4 + if obcndx==1 + suffix='_south'; % SOUTH + Dstp=squeeze(h(1,:)); + hrm=2*squeeze(Awave(l,1,:)); + cdir=3*pi/2-deg2rad(squeeze(Dwave(l,1,:)))-grid_angle; + cfrq=2*pi./squeeze(Pwave(l,1,:)); + elseif obcndx==2 + suffix='_east'; % EAST + Dstp=squeeze(h(:,L)); + hrm=2*squeeze(Awave(l,:,L)); + cdir=3*pi/2-deg2rad(squeeze(Dwave(l,:,L)))-grid_angle; + cfrq=2*pi./squeeze(Pwave(l,:,L)); + elseif obcndx==3 + suffix='_north'; % NORTH + Dstp=squeeze(h(M,:)); + hrm=2*squeeze(Awave(l,M,:)); + cdir=3*pi/2-deg2rad(squeeze(Dwave(l,M,:)))-grid_angle; + cfrq=2*pi./squeeze(Pwave(l,M,:)); + elseif obcndx==4 + suffix='_west'; % WEST + Dstp=squeeze(h(:,1)); + hrm=2*squeeze(Awave(l,:,1)); + cdir=3*pi/2-deg2rad(squeeze(Dwave(l,:,1)))-grid_angle; + cfrq=2*pi./squeeze(Pwave(l,:,1)); + end + dd = Dstp'; %+Tide(l) + khd = dd.*cfrq.*cfrq./g; + kh = sqrt( khd.*khd + khd./(1.0 + khd.*(0.6666666666 ... + +khd.*(0.3555555555 + khd.*(0.1608465608 ... + +khd.*(0.0632098765 + khd.*(0.0217540484 ... + +khd.*0.0065407983)))))) ); + kw=kh./dd; + cosw=cos(cdir); + sinw=sin(cdir); + wac=0.125*g*(hrm.^2)./cfrq; + wkx=kw.*cosw; + wke=kw.*sinw; + nc{['wac' suffix]}(l,:)=wac; + nc{['wkx' suffix]}(l,:)=wkx; + nc{['wke' suffix]}(l,:)=wke; + end + end + close(nc); + + end % M +end % Y +% +% End +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Coupling_tools/CROCO/prepro_soda.m b/Coupling_tools/CROCO/prepro_soda.m index 3ba9797bd55287e75ffa411903ae60496b3cf6e5..2bd90cae5a0179c9165783b3897079ed8cffaa67 100644 --- a/Coupling_tools/CROCO/prepro_soda.m +++ b/Coupling_tools/CROCO/prepro_soda.m @@ -1,4 +1,4 @@ clear all close all start ; -make_OGCM ; +make_OGCM_SODA ; diff --git a/Forecast_tools/download_ECCO_frcst.m b/Forecast_tools/download_ECCO_frcst.m index 64653c233bad651bf8f3ff98ec828f823cab2db6..c0f3703ce4c978316ab46bf1ab8a2c052b7c2eae 100644 --- a/Forecast_tools/download_ECCO_frcst.m +++ b/Forecast_tools/download_ECCO_frcst.m @@ -161,13 +161,13 @@ if ~exist(ecco_name) % % Extract ECCO % - extract_ECCO_frcst(FRCST_dir,FRCST_prefix,prefix,suffix,tindex,missval,... + extract_ECCO_frcst(FRCST_dir,FRCST_prefix,prefix,suffix,tindex,missval,... lon,lon_u,lat,lat_v,depth,... krange,jrange,jrange_v,... i1min,i1max,i2min,i2max,i3min,i3max,... i1min_u,i1max_u,i2min_u,i2max_u,i3min_u,i3max_u,... time,Yorig) else - disp(' ECCO file allready exist') + disp(' ECCO file already exists') end return diff --git a/Forecast_tools/download_GFS.m b/Forecast_tools/download_GFS.m index 6060cd28ac5acc6a170772898030546b54407f75..ec74fdc261a3279fe19b8b6eedc1d695544fa6ca 100644 --- a/Forecast_tools/download_GFS.m +++ b/Forecast_tools/download_GFS.m @@ -193,7 +193,8 @@ for ihind=1:4*hdays % loop on files until time=yesterday 12Z radlw(n,:,:),radlw_in(n,:,:),radsw(n,:,:)]=... get_GFS(fname,mask,tndxdap,jrange,i1min,i1max,i2min,i2max, ... i3min,i3max,missvalue); - end + end + %pause(20) end % ihind (nb of hindcast record) % %================================================================== @@ -255,6 +256,7 @@ for tndx=t1+1:it:tend radlw(n,:,:),radlw_in(n,:,:),radsw(n,:,:)]=... get_GFS(fname,mask,tndxdap,jrange,i1min,i1max,i2min,i2max, ... i3min,i3max,missvalue); + %pause(20) end % %================================================================== diff --git a/Forecast_tools/download_mercator_frcst_python.m b/Forecast_tools/download_mercator_frcst_python.m index 9545df18c7a093e0cfff28faad98ba78916d6145..aa11005df5d5076375c7b3f4b10953ed796a949a 100644 --- a/Forecast_tools/download_mercator_frcst_python.m +++ b/Forecast_tools/download_mercator_frcst_python.m @@ -65,20 +65,18 @@ end rundate_str=date; date_run=datenum(rundate_str); rundate=date_run-datenum(Yorig,1,1); -% -% We consider for hdays=2 (hindcast NRT analysis + hindcast simulation P0) for i=1:lh+1 - time1(i)=date_run-(lh-i+1); + time1(i)=datenum(rundate_str)-(lh+2-i); end -time2=date_run; -for j=1:lf - time3(j)=date_run+j-1; % because P1 is date_run +time2=datenum(rundate_str); +for j=1:lf+1 + time3(j)=datenum(rundate_str)+j; end time=cat(2,time1,time2,time3); tiempo_inicial = time1(1); +tiempo_final = time3(end); tiempo_p0 = time1(lh); tiempo_p1 = date_run; -tiempo_final = time3(end); % if (lonmin > 180) lonmin = lonmin - 360; @@ -97,12 +95,12 @@ disp([' ']) if download_raw_data % disp(['--> Extraction Forecasts from operational date run (now) : ', datestr(tiempo_p1,'yyyy-mm-dd')]) - if lh ~= 1 + if lh ~= 0 for i=1:lh-1 - disp(['Get Hindcasts NRT ana.: ',datestr(time1(i),'yyyy-mm-dd')]) + disp(['Get Hindcasts NRT ana.: ',datestr(time1(i),'yyyy-mm-dd')]) end - end - disp(['Get Hindcasts P0 : ',datestr(tiempo_p0,'yyyy-mm-dd')]) + disp(['Get Hindcasts P0 : ',datestr(tiempo_p0,'yyyy-mm-dd')]) + end disp(['Get ',num2str(lf),' Forecasts : from ',datestr(tiempo_p1,'yyyy-mm-dd'),' to ',datestr(tiempo_final, 'yyyy-mm-dd')]) % % diff --git a/Forecast_tools/get_GFS.m b/Forecast_tools/get_GFS.m index 667463e388d902d05f92cffe7288e8e268568684..628590b64df5a2624f5778557d2a3fa05c480a23 100644 --- a/Forecast_tools/get_GFS.m +++ b/Forecast_tools/get_GFS.m @@ -47,29 +47,32 @@ t=readdap(fname,'time',trange); disp(['TRANGE=',num2str(trange)]) disp(['GFS raw time=',sprintf('%5.3f',t)]) t=t+365; % put it in "matlab" time -disp(['GFS: ',datestr(t)]) +disp(['GFS: ',datestr(datenum(t),'yyyy-mm-dd HH:MM:SS')]) disp('====================================================') %disp('u...') u=mask.*getdap('',fname,'ugrd10m',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); u(abs(u)>=missvalue)=NaN; +pause(1) %disp('v...') v=mask.*getdap('',fname,'vgrd10m',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); v(abs(v)>=missvalue)=NaN; +pause(1) %disp('ty...') ty=mask.*getdap('',fname,'vflxsfc',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); ty(abs(ty)>=missvalue)=NaN; +pause(1) %disp('tx...') tx=mask.*getdap('',fname,'uflxsfc',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); tx(abs(tx)>=missvalue)=NaN; - +pause(1) %disp('skt...') %skt=mask.*getdap('',fname,'tmpsfc',trange,'',jrange,... @@ -80,36 +83,43 @@ tx(abs(tx)>=missvalue)=NaN; tair=mask.*getdap('',fname,'tmp2m',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); tair(abs(tair)>=missvalue)=NaN; +pause(1) %disp('rhum...') rhum=mask.*getdap('',fname,'rh2m',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); rhum(abs(rhum)>=missvalue)=NaN; +pause(1) %disp('prate...') prate=mask.*getdap('',fname,'pratesfc',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); prate(abs(prate)>=missvalue)=NaN; +pause(1) %disp('down radlw...') dradlw=mask.*getdap('',fname,'dlwrfsfc',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); dradlw(abs(dradlw)>=missvalue)=NaN; +pause(1) %disp('up radlw...') uradlw=mask.*getdap('',fname,'ulwrfsfc',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); uradlw(abs(uradlw)>=missvalue)=NaN; +pause(1) %disp('down radsw...') dradsw=mask.*getdap('',fname,'dswrfsfc',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); dradsw(abs(dradsw)>=missvalue)=NaN; +pause(1) %disp('up radsw...') uradsw=mask.*getdap('',fname,'uswrfsfc',trange,'',jrange,... i1min,i1max,i2min,i2max,i3min,i3max); uradsw(abs(uradsw)>=missvalue)=NaN; +pause(1) % % Transform the variables % diff --git a/Forecast_tools/make_GFS.m b/Forecast_tools/make_GFS.m index 967a7048d566a03c968f475e2b44939ce90ab740..4d956884a3aac3b646640d69c9340023ecc3c159 100644 --- a/Forecast_tools/make_GFS.m +++ b/Forecast_tools/make_GFS.m @@ -95,6 +95,10 @@ if Download_data==1 % Download data with DODS (the download matlab routine depends on the OGCM) % disp('Download data...') + disp('--> NB: if you receive this message during download : ') + disp(' " syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR ..."') + disp(' " you maybe reach the rate limit authorized, there is a ban on your IP due to an “Over Rate Limit†violation, try to wait more than 10 minutes before to do it again') + disp('') download_GFS(today,lonmin,lonmax,latmin,latmax,FRCST_dir,Yorig,it) % end diff --git a/Forecast_tools/make_OGCM_frcst.m b/Forecast_tools/make_OGCM_frcst.m index fabf675623355fbf006bd075b041808f50e363a9..e6400335944aaba6f57ec73ea7735eab4b17176a 100644 --- a/Forecast_tools/make_OGCM_frcst.m +++ b/Forecast_tools/make_OGCM_frcst.m @@ -56,37 +56,27 @@ if strcmp(OGCM,'mercator') for i=1:hdays+1 time1(i)=date_run-(hdays+1-i); end - rundate=datestr(time1(1+hdays),'yyyy-mm-ddT00:00:00Z'); + date_deb=datestr(time1(1+hdays),'yyyy-mm-ddT00:00:00Z'); if hdays ~= 0 - % date_deb = rundate - hindcast days : example 2024-02-28T00:00:00Z-PT2D - % for a validity date 2 days before run 20240228 00hUTC ... + % date_deb = rundate - hindcast days : example 2024-02-28T00:00:00Z-PT1D + % for a validity date 1 day before run 20240228 00hUTC ... date_deb=strcat(datestr(time1(1+hdays),'yyyy-mm-ddT00:00:00') ,['-PT',num2str(hdays),'D']); - else - date_deb=rundate; end end % % Set generic OGCM file name % FRCST_prefix=[OGCM,'_']; -OGCM_name=[FRCST_dir,FRCST_prefix,num2str(rundate),'.cdf']; +OGCM_name=[FRCST_dir,FRCST_prefix,num2str(date_deb),'.cdf']; % -if strcmp(OGCM,'ECCO') - % - % ECCO DODS URL - % - % Kalman filter - % - url = 'http://ecco.jpl.nasa.gov/thredds/dodsC/las/kf080/kf080_'; - % -elseif strcmp(OGCM,'mercator') +if strcmp(OGCM,'mercator') % % MERCATOR : see get_file_python_mercator.m % - raw_mercator_name=[FRCST_dir,'raw_mercator_',num2str(rundate),'.nc']; - mercator_name=[FRCST_dir,FRCST_prefix,num2str(rundate),'.cdf']; + raw_mercator_name=[FRCST_dir,'raw_mercator_',num2str(date_deb),'.nc']; + mercator_name=[FRCST_dir,FRCST_prefix,num2str(date_deb),'.cdf']; else - error(['Unknown OGCM: ',OGCM]) + error(['Unknown OGCM for forecast: ',OGCM]) end % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -125,13 +115,7 @@ if Download_data % % Download data (matlab routine for download depends on OGCM) % - if strcmp(OGCM,'ECCO') - disp('Download data...') - eval(['OGCM_name=download_', ... - OGCM,'_frcst(lonmin,lonmax,latmin,latmax,',... - 'FRCST_dir,FRCST_prefix,url,Yorig);']) - - elseif strcmp(OGCM,'mercator') + if strcmp(OGCM,'mercator') % % Use of Copernicusmarine client (python) % -> script download_mercator_frcst_python.m @@ -143,9 +127,11 @@ if Download_data 'hdays,fdays,', ... 'lonmin,lonmax,latmin,latmax,hmax,', ... 'FRCST_dir,FRCST_prefix,mercator_name,Yorig);']) + + else + error(['Unknown OGCM for forecast: ',OGCM]) end end - %--------------------------------------------------------------- % Get OGCM grid %--------------------------------------------------------------- @@ -160,15 +146,10 @@ Z=-nc{'depth'}(:); NZ=length(Z); NZ=NZ-rmdepth; Z=Z(1:NZ); - %--------------------------------------------------------------- % Get time array %--------------------------------------------------------------- -if strcmp(OGCM,'ECCO') - time=[90 270]; - time_cycle=360; - trange=[1 1]; -elseif strcmp(OGCM,'mercator') +if strcmp(OGCM,'mercator') OGCM_time=nc{'time'}(:); time_cycle=0; delta=1; % >1 if subsampling @@ -177,15 +158,16 @@ elseif strcmp(OGCM,'mercator') for i=1:length(trange) time(i)=OGCM_time(trange(i)); end +else + error(['Unknown OGCM for forecast: ',OGCM]) end close(nc) - %--------------------------------------------------------------- % Initial file %--------------------------------------------------------------- if makeini==1 % - ininame=[ini_prefix,num2str(date_deb),nc_suffix]; + ininame=[ini_prefix,num2str(rundate),nc_suffix]; disp(['-> Create an initial file for date ',num2str(date_deb);]) create_inifile(ininame,grdname,CROCO_title,... theta_s,theta_b,hc,N,... @@ -199,19 +181,18 @@ if makeini==1 close(nc_ini) if hdays ~= 0 if (isoctave == 0) - eval(['!cp ',ininame,' ',ini_prefix,'hct_',num2str(date_deb),nc_suffix]) + eval(['!cp ',ininame,' ',ini_prefix,'hct_',num2str(rundate),nc_suffix]) else - eval(['cp ',ininame,' ',ini_prefix,'hct_',num2str(date_deb),nc_suffix]) + eval(['cp ',ininame,' ',ini_prefix,'hct_',num2str(rundate),nc_suffix]) end end end - %--------------------------------------------------------------- % Clim and Bry files %--------------------------------------------------------------- if makeclim==1 | makebry==1 if makebry==1 - bryname=[bry_prefix,num2str(date_deb),nc_suffix]; + bryname=[bry_prefix,num2str(rundate),nc_suffix]; disp([' ']) disp(['-> Create bry files from ',num2str(date_deb);]) create_bryfile(bryname,grdname,CROCO_title,[1 1 1 1],... @@ -222,7 +203,7 @@ if makeclim==1 | makebry==1 nc_bry=[]; end if makeclim==1 - clmname=[clm_prefix,num2str(date_deb),nc_suffix]; + clmname=[clm_prefix,num2str(rundate),nc_suffix]; disp(['-> Create clim files from ',num2str(date_deb);]) create_climfile(clmname,grdname,CROCO_title,... theta_s,theta_b,hc,N,... @@ -231,7 +212,6 @@ if makeclim==1 | makebry==1 else nc_clm=[]; end - %--------------------------------------------------------------- % Perform interpolations for all selected records %--------------------------------------------------------------- @@ -242,7 +222,6 @@ for tndx=1:length(time) nc_clm,nc_bry,lon,lat,angle,h,pm,pn,rmask, ... tndx,vtransform,obc) end - % % Close CROCO files % @@ -261,21 +240,21 @@ if makeplot==1 disp(' ') disp(' Make a few plots...') if makeini==1 - ininame=[ini_prefix,num2str(date_deb),nc_suffix]; + ininame=[ini_prefix,num2str(rundate),nc_suffix]; figure test_clim(ininame,grdname,'temp',1,coastfileplot) figure test_clim(ininame,grdname,'salt',1,coastfileplot) end if makeclim==1 - clmname=[clm_prefix,num2str(date_deb),nc_suffix]; + clmname=[clm_prefix,num2str(rundate),nc_suffix]; figure test_clim(clmname,grdname,'temp',1,coastfileplot) figure test_clim(clmname,grdname,'salt',1,coastfileplot) end if makebry==1 - bryname=[bry_prefix,num2str(date_deb),nc_suffix]; + bryname=[bry_prefix,num2str(rundate),nc_suffix]; figure test_bry(bryname,grdname,'temp',1,obc) figure diff --git a/Forecast_tools/plot_forecast_croco.m b/Forecast_tools/plot_forecast_croco.m index 82375860ceed13c236559114fe34acbdbb54aeb9..ffaad1614135c6257ace056ed17202d6f332c1ec 100644 --- a/Forecast_tools/plot_forecast_croco.m +++ b/Forecast_tools/plot_forecast_croco.m @@ -16,6 +16,7 @@ dirout = './Figs/'; eval(['!mkdir ',dirout]) grdname = [dirin,'SCRATCH/croco_grd.nc']; frcname = [dirin,'SCRATCH/croco_frc_GFS_0.nc']; +hisname = [dirin,'SCRATCH/croco_his.nc']; avgname = [dirin,'SCRATCH/croco_avg.nc']; coastfile = [dirin,'coastline_l.mat']; @@ -26,7 +27,7 @@ skip_cur = 2; zoom = 0; fsize = 14; % Font size nx = 3; % number of days -titlemark = ['IRD: ',datestr(now)]; +titlemark = ['CROCO: ',datestr(now)]; mean_currents = 0; % 1: plot 0-500 mean currents % 0: surface currents @@ -36,6 +37,8 @@ plot_fig_2 = 0; % surface currents plot_fig_3 = 1; % SST + surface currents plot_fig_4 = 1; % SSH plot_fig_5 = 0; % HBL +plot_fig_6 = 0; % max SSH (his file) + % offset to convert UTC model time % into local matlab time (days] @@ -43,7 +46,7 @@ days_off = timezone/24+datenum(Yorig,1,1); % %************ End of users's defined parameter ************ % -% time (in matlab time) +/bin/bash: q : commande introuvable % today=floor(now); % @@ -75,7 +78,6 @@ dzv=zv(2:end,:,:)-zv(1:end-1,:,:); % coastline %bounds=[lonmin lonmax latmin latmax]; res='l'; %coastfile=make_coast2(bounds,res) - % barwidth=0.5; barheight=0.04; @@ -437,3 +439,55 @@ for tndx=hdays+1:nx+hdays end end end % plot_fig_5 + +if plot_fig_6, +% +% 6: Sea Surface Height max +% + disp('Plot SSH max...') + itime=0; + cmin=0; cint=0.5; cmax=7; + + outname=['MaxSSH']; + figure('position',[5 5 1000 600]); + subplot('position',[0 0.2 1 .75]); + + nc=netcdf(hisname); + N=length(nc('s_rho')); + time=nc{'scrum_time'}(:)/86400+days_off; + t1=datestr(time(1),1); + t2=datestr(time(end),1); + ssh=squeeze(max(nc{'zeta'}(:))); + mask=nc{'mask_rho'}(:); + ssh=ssh.*mask./mask; + close(nc) + ssh=get_missing_val(lon,lat,ssh); + + m_proj('mercator','lon',[lonmin lonmax],'lat',[latmin latmax]); + [C0,h0]=m_contourf(lon,lat,ssh,[cmin:cint:cmax]); + clabel(C0,h0,'fontsize',14,'LabelSpacing',100) + caxis([cmin cmax]) + m_usercoast(coastfile,'patch',[.0 .0 .0]); + m_grid('box','fancy','tickdir','in','FontSize',fsize-2); + set(findobj('tag','m_grid_color'),'facecolor','white') + title(['Max Elevation [',t1,' - ',t2,']'],'FontSize',fsize+3) + ht=m_text(lonmax-0.05,latmax-.05,titlemark); + set(ht,'horizontalalignment','right','fontsize',fsize-3,'Color','w'); + set(gca,'Layer','top'); + + subplot('Position',[0.5-0.5*barwidth vmargin barwidth barheight]) + x=[0:1]; + y=[cmin:cint:cmax]; + [X,Y]=meshgrid(x,y); + caxis([cmin cmax]) + contourf(Y,X,Y,y) + set(gca,'XTick',[cmin:cint:cmax],'YTickLabel',[' ']) + set(gca,'FontSize',fsize) + set(gca,'Layer','top'); + xlabel('SSH [m]','FontSize',fsize) + set(gcf, 'PaperPositionMode', 'auto'); + + export_fig -transparent file.pdf + eval(['! mv -f file.pdf ',dirout,outname,'.pdf']); +end % plot_fig_6 + diff --git a/Oforc_OGCM/Copernicus_Marine_Toolbox_installation.md b/Oforc_OGCM/Copernicus_Marine_Toolbox_installation.md index 2820c4b2140b7ce19abd6ee8ce52aa6fa51e9f41..685502c405c9f7aa6dc7c9bceec8a1fe3e1c1076 100644 --- a/Oforc_OGCM/Copernicus_Marine_Toolbox_installation.md +++ b/Oforc_OGCM/Copernicus_Marine_Toolbox_installation.md @@ -16,13 +16,21 @@ Prerequisites: Note that you can use micromamba instead of mamba to install the python environment you can download also the file copernicusmarine_env.yml (See Instructions on 2) - -To use it, you need to activate the cmt_1.0 : +Firstly, you need to activate the environment cmt_1.0 : > mamba activate cmt_1.0 -You will have access to the copernicusmarine executable with various sub-command, very useful to get, extract and download data from the Copernicus Marine Data Store. +The location of the executable (here after pathCMC) will be found typing : + > ls $CONDA_PREFIX/bin/copernicusmarine + + Note that the value returned by your terminal here, will be your pathCMC to fill in your crocotools_param.m + and download_glorys_data.sh script + +You will have access to the copernicusmarine executable with various sub-command, very useful to get, extract and download data from the Copernicus Marine Data Store : > copernicusmarine -h +Or alternatively: + > $CONDA_PREFIX/bin/copernicusmarine -h + To be used in the Matlab croco_tools, in the crocotools_param.m, you need to : diff --git a/Oforc_OGCM/create_OGCM.m b/Oforc_OGCM/create_OGCM.m index 24e4ae1e8f4ec52bcc54ac3e318122b3fbacfa05..88f0469002ac1d2a37c6e9b7521ee2bb1d33b7b0 100644 --- a/Oforc_OGCM/create_OGCM.m +++ b/Oforc_OGCM/create_OGCM.m @@ -80,6 +80,18 @@ nc{'vbar'}.long_name='MERIDIONAL BAROTROPIC VELOCITY'; nc{'vbar'}.units=ncchar('m/sec'); nc{'vbar'}.units='m/sec'; nc{'vbar'}.missing_value=missval; +%nc{'taux'}=ncfloat('time','latU','lonU') ; +%nc{'taux'}.long_name=ncchar('TAU_X'); +%nc{'taux'}.long_name='TAU_X'; +%nc{'taux'}.units=ncchar('N.m-2'); +%nc{'taux'}.units='N.m-2'; +%nc{'taux'}.missing_value=missval; +%nc{'tauy'}=ncfloat('time','latV','lonV') ; +%nc{'tauy'}.long_name=ncchar('TAU_Y'); +%nc{'tauy'}.long_name='TAU_Y'; +%nc{'tauy'}.units=ncchar('N.m-2'); +%nc{'tauy'}.units='N.m-2'; +%nc{'tauy'}.missing_value=missval; nc{'ssh'}=ncfloat('time','latT','lonT') ; nc{'ssh'}.long_name=ncchar('SEA LEVEL HEIGHT'); nc{'ssh'}.long_name='SEA LEVEL HEIGHT'; @@ -129,6 +141,8 @@ nc{'time'}(tndx)=time(tndx); % if length(time)==1 nc{'ssh'}(tndx,:,:)=ssh; +% nc{'taux'}(tndx,:,:)=taux; +% nc{'tauy'}(tndx,:,:)=tauy; u1=u; v1=v; nc{'u'}(tndx,:,:,:)=u1; @@ -137,6 +151,8 @@ nc{'time'}(tndx)=time(tndx); nc{'salt'}(tndx,:,:,:)=salt; else nc{'ssh'}(tndx,:,:)=squeeze(ssh(tndx,:,:)); +% nc{'taux'}(tndx,:,:)=squeeze(taux(tndx,:,:)); +% nc{'tauy'}(tndx,:,:)=squeeze(tauy(tndx,:,:)); u1=squeeze(u(tndx,:,:,:)); v1=squeeze(v(tndx,:,:,:)); nc{'u'}(tndx,:,:,:)=u1; diff --git a/Oforc_OGCM/download_glorys_data.sh b/Oforc_OGCM/download_glorys_data.sh index 819b4bb7d2b0bcf2f2cd6e4ed405a5e168e93dad..8b20c40a142649eef765cf9a5ea86858ecf95afb 100755 --- a/Oforc_OGCM/download_glorys_data.sh +++ b/Oforc_OGCM/download_glorys_data.sh @@ -1,6 +1,6 @@ #!/bin/bash set -e -source ../../myenv_mypath.sh +#source ../../myenv_mypath.sh ########################## # python # ====== @@ -15,9 +15,9 @@ source ../../myenv_mypath.sh user='XXXXXXXXX' password='XXXXXXXXX' ### -### motu client ### -path_to_motu="${OCE}/../../croco_tools" # Use croco's motuclient croco_tools/Forecast_tools -# if you want to use your own motu, leave empty +### Copernicus marine client ### +# this is the pathCMC (described on the readme file Copernicus_Marine_Toolbox_installation.md) +path_to_cmt="/home/to/your/bin/copernicusmarine" ### YEAR_START=2017 MONTH_START=1 @@ -36,7 +36,8 @@ DT_TIME= # Value in days. Default value when empty is 1 month. kdata="MONTHLY" # DAILY or MONTHLY #################### READ_GRD=1 # read croco_grid to find lon/lat min-max -INPUT_GRD="${OCE_FILES_DIR}/croco_grd.nc" +#INPUT_GRD="${OCE_FILES_DIR}/croco_grd.nc" +INPUT_GRD="CROCO_FILES/croco_grd.nc" # In case you don't want to read croco_grd # Please use -180/180 -90/90 format lon_min="-90" @@ -59,16 +60,14 @@ for field in ${vars}; do done ### -### motu command ### -if [[ -z ${path_to_motu} ]]; then - command_line="python -m motuclient" +if [[ -z ${path_to_cmt} ]]; then + echo "You need to define your path_to_cmt variable ..." + exit 12 else - command_line="${path_to_motu}/Forecast_tools/motuclient-python/motuclient.py" + command_line=${path_to_cmt}" " fi -## motu info -motu_url_reana='http://my.cmems-du.eu/motu-web/Motu' -service_id_reana='GLOBAL_MULTIYEAR_PHY_001_030-TDS' +## cmt info if [[ ${kdata} == "DAILY" ]]; then product_id_reana='cmems_mod_glo_phy_my_0.083_P1D-m' elif [[ ${kdata} == "MONTHLY" ]]; then @@ -113,18 +112,20 @@ for YEAR in `seq ${YEAR_START} ${YEAR_END}`; do daystrt=1 start_date=$( printf "%04d-%02d-%02d" $YEAR $MONTH $daystrt) end_date=$(printf `date +"%Y-%m-%d" -d "${start_date} +1 month - 1 day"`) - OUTNAME="${OUTDIR}/raw_motu_${PREFIX}_Y${YEAR}M${MONTH}.nc" - ${command_line} --motu ${motu_url_reana} --service-id ${service_id_reana} --product-id ${product_id_reana} --longitude-min ${lon_min} --longitude-max ${lon_max} --latitude-min ${lat_min} --latitude-max ${lat_max} --date-min "${start_date} 12:00:00" --date-max "${end_date} 12:00:00" --depth-min 0.493 --depth-max 5727.918 ${variables} --out-dir ./ --out-name ${OUTNAME} --user ${user} --pwd ${password} + OUTNAME="${OUTDIR}/raw_${PREFIX}_Y${YEAR}M${MONTH}.nc" + echo "${command_line} subset -i ${product_id_reana} -x ${lon_min} -X ${lon_max} -y ${lat_min} -Y ${lat_max} -t "${start_date}" -T "${end_date}" -z 0.493 -Z 5727.918 ${variables} -o ./ -f ${OUTNAME} --username ${user} --password ${password} --force-download" + ${command_line} subset -i ${product_id_reana} -x ${lon_min} -X ${lon_max} -y ${lat_min} -Y ${lat_max} -t "${start_date}" -T "${end_date}" -z 0.493 -Z 5727.918 ${variables} -o ./ -f ${OUTNAME} --username ${user} --password ${password} --force-download elif [[ ${kdata} == "DAILY" ]]; then if [ -z ${DT_TIME} ]; then start_date=$( printf "%04d-%02d-%02d" $YEAR $MONTH $daystrt) end_date=$(printf `date +"%Y-%m-%d" -d "${start_date} +1 month - 1 day"`) - OUTNAME="${OUTDIR}/raw_motu_${PREFIX}_Y${YEAR}M${MONTH}.nc" - ${command_line} --motu ${motu_url_reana} --service-id ${service_id_reana} --product-id ${product_id_reana} --longitude-min ${lon_min} --longitude-max ${lon_max} --latitude-min ${lat_min} --latitude-max ${lat_max} --date-min "${start_date} 12:00:00" --date-max "${end_date} 12:00:00" --depth-min 0.493 --depth-max 5727.918 ${variables} --out-dir ./ --out-name ${OUTNAME} --user ${user} --pwd ${password} + OUTNAME="${OUTDIR}/raw_${PREFIX}_Y${YEAR}M${MONTH}.nc" + echo "${command_line} subset -i ${product_id_reana} -x ${lon_min} -X ${lon_max} -y ${lat_min} -Y ${lat_max} -t "${start_date}" -T "${end_date}" -z 0.493 -Z 5727.918 ${variables} -o ./ -f ${OUTNAME} --username ${user} --password ${password} --force-download" + ${command_line} subset -i ${product_id_reana} -x ${lon_min} -X ${lon_max} -y ${lat_min} -Y ${lat_max} -t "${start_date}" -T "${end_date}" -z 0.493 -Z 5727.918 ${variables} -o ./ -f ${OUTNAME} --username ${user} --password ${password} --force-download else for DAY in `seq ${dstart} ${DT_TIME} ${dend}`; do - tmpoutname="${OUTDIR}/raw_motu_${PREFIX}_Y${YEAR}M${MONTH}D${DAY}.nc" + tmpoutname="${OUTDIR}/raw_${PREFIX}_Y${YEAR}M${MONTH}D${DAY}.nc" start_date=$( printf "%04d-%02d-%02d" $YEAR $MONTH $DAY ) end_date=$( date +"%Y-%m-%d" -d "${start_date} +${DT_TIME} day - 1 day" ) tmpsdate=$( echo `date -d ${start_date} +"%Y%m"` ) @@ -137,12 +138,13 @@ for YEAR in `seq ${YEAR_START} ${YEAR_END}`; do nbdays=$(( (${second_date} - ${first_date})/86400 )) end_date=$( date +"%Y-%m-%d" -d "${end_date} - ${nbdays} day" ) fi - ${command_line} --motu ${motu_url_reana} --service-id ${service_id_reana} --product-id ${product_id_reana} --longitude-min ${lon_min} --longitude-max ${lon_max} --latitude-min ${lat_min} --latitude-max ${lat_max} --date-min "${start_date} 12:00:00" --date-max "${end_date} 12:00:00" --depth-min 0.493 --depth-max 5727.918 ${variables} --out-dir ./ --out-name ${tmpoutname} --user ${user} --pwd ${password} + echo "${command_line} subset -i ${product_id_reana} -x ${lon_min} -X ${lon_max} -y ${lat_min} -Y ${lat_max} -t "${start_date}" -T "${end_date}" -z 0.493 -Z 5727.918 ${variables} -o ./ -f ${tmpoutname} --username ${user} --password ${password} --force-download" + ${command_line} subset -i ${product_id_reana} -x ${lon_min} -X ${lon_max} -y ${lat_min} -Y ${lat_max} -t "${start_date}" -T "${end_date}" -z 0.493 -Z 5727.918 ${variables} -o ./ -f ${tmpoutname} --username ${user} --password ${password} --force-download ncks -O -F --mk_rec_dmn time ${tmpoutname} ${tmpoutname} done cd ${OUTDIR} - ncrcat -O raw_motu_${PREFIX}_Y${YEAR}M${MONTH}D*.nc raw_motu_${PREFIX}_Y${YEAR}M${MONTH}.nc - rm -r raw_motu_${PREFIX}_Y${YEAR}M${MONTH}D*.nc + ncrcat -O raw_${PREFIX}_Y${YEAR}M${MONTH}D*.nc raw_${PREFIX}_Y${YEAR}M${MONTH}.nc + rm -r raw_${PREFIX}_Y${YEAR}M${MONTH}D*.nc cd - fi fi diff --git a/Oforc_OGCM/extract_SODA.m b/Oforc_OGCM/extract_SODA.m index d9ac1d4f627e6aa08f9425808e30e238ff5b28e3..1790f396155e64cf27ac055e74bad01c25db2f0e 100644 --- a/Oforc_OGCM/extract_SODA.m +++ b/Oforc_OGCM/extract_SODA.m @@ -152,6 +152,6 @@ end % temp,salt,u,v,ssh,taux,tauy,Yorig) create_OGCM([SODA_dir,SODA_prefix,'Y',num2str(year),'M',num2str(month),'.cdf'],... lon,lat,lon,lat,lon,lat,depth,time,... - temp,salt,u,v,ssh,Yorig) + temp,salt,u,v,ssh,Yorig) % return diff --git a/Oforc_OGCM/make_OGCM_SODA.m b/Oforc_OGCM/make_OGCM_SODA.m index b653a944444b33b9edb50689b676a560b2c6556d..f59e5e54c0f8ec650b2a2044832c2bf95b552635 100644 --- a/Oforc_OGCM/make_OGCM_SODA.m +++ b/Oforc_OGCM/make_OGCM_SODA.m @@ -300,10 +300,10 @@ if makeclim==1 | makebry==1 Yp=Y; for aa=1:itolap_p tndx_OGCM(aa)=ntimes; - end + end else for aa=1:itolap_p - tndx_OGCM(aa)=aa; + tndx_OGCM(aa)=aa; end; end % diff --git a/Oforc_OGCM/prepro_soda.m b/Oforc_OGCM/prepro_soda.m index 3ba9797bd55287e75ffa411903ae60496b3cf6e5..2bd90cae5a0179c9165783b3897079ed8cffaa67 100644 --- a/Oforc_OGCM/prepro_soda.m +++ b/Oforc_OGCM/prepro_soda.m @@ -1,4 +1,4 @@ clear all close all start ; -make_OGCM ; +make_OGCM_SODA ; diff --git a/Oforc_OGCM/write_mercator_multi.m b/Oforc_OGCM/write_mercator_multi.m index 717629e9c719b0587086692ca39798f0c24c8991..3f4878da038d8f5e4efdca245357782603e3a31f 100644 --- a/Oforc_OGCM/write_mercator_multi.m +++ b/Oforc_OGCM/write_mercator_multi.m @@ -51,7 +51,6 @@ lat = nc{'latitude'}(:); depth = nc{'depth'}(:); time = nc{'time'}(:); time = time / 24 + datenum(1950,1,1) - datenum(Yorig,1,1); - % % Get SSH % @@ -139,7 +138,6 @@ close(nc) % close raw_mercator_name % % Create the Mercator file % - create_OGCM([OGCM_dir,OGCM_prefix,thedatemonth,'.cdf'],... lon,lat,lon,lat,lon,lat,depth,time,... squeeze(temp),squeeze(salt),squeeze(u),... diff --git a/Preprocessing_tools/make_config.m b/Preprocessing_tools/make_config.m index d6e53acfd725406bd8d52759e8a505097c1f2ef4..94b89ab97f69c3e85168aa23229b533262110466 100644 --- a/Preprocessing_tools/make_config.m +++ b/Preprocessing_tools/make_config.m @@ -47,5 +47,5 @@ make_clim % Interannual forcing % %make_NCEP -%make_OGCM +%make_OGCM_mercator (SODA) return diff --git a/crocotools_param.m b/crocotools_param.m index e765243b75d2977bfdc6f7269032b95b29ba741d..28bf7d640e272d02868a183b48a16eb754526e82 100644 --- a/crocotools_param.m +++ b/crocotools_param.m @@ -5,7 +5,7 @@ % % This file is used by make_grid.m, make_forcing.m, % make_clim.m, make_biol.m, make_bry.m, make_tides.m, -% make_NCEP.m, make_OGCM.m, make_... +% make_NCEP.m, make_OGCM_*.m, make_... % % Further Information: % http://www.croco-ocean.org @@ -221,7 +221,7 @@ pathfinder_sst_name=[DATADIR,... % % 4 - Open boundaries and initial conditions parameters % used by make_clim.m, make_biol.m, make_bry.m -% make_OGCM.m and make_OGCM_frcst.m +% make_OGCM_*.m and make_OGCM_frcst.m % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % @@ -235,7 +235,7 @@ zref = -1000; % % initial/boundary data options (1 = process) % (used in make_clim, make_biol, make_bry, -% make_OGCM.m and make_OGCM_frcst.m) +% make_OGCM_*.m and make_OGCM_frcst.m) % makeini = 1; % initial data makeclim = 1; % climatological data (for boundaries and nudging layers) @@ -322,7 +322,7 @@ Z0 = 1; % Mean depth of tide gauge %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % 6 - Reference date and simulation times -% (used for make_tides, make_CFSR (or make_NCEP), make_OGCM) +% (used for make_tides, make_CFSR (or make_NCEP), make_OGCM_*) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % @@ -401,9 +401,9 @@ QSCAT_clim_file = [DATADIR,'QuikSCAT_clim/',... % QSCAT climatolog % Reformat_ECMWF = 1; ECMWF_dir = [FORC_DATA_DIR,'ECMWF_',CROCO_config,'/']; % ERA-I data dir. [croco format] -My_ECMWF_dir = [FORC_DATA_DIR,'ERAI/']; % ERA-I native data downloaded - % with python script -itolap_ecmwf = 3; % 3 records for daily ECMWF +My_ECMWF_dir = [FORC_DATA_DIR,'ERAI/']; % ERA-I native data downloaded + % with python script +itolap_ecmwf = 3; % 3 records for daily ECMWF % %-------------------------------------------------- % Options for make_ERA5 @@ -417,7 +417,7 @@ itolap_era5 = 2; % 2 records = 2 h % % %-------------------------------------------- -% Options for make_OGCM or make_OGCM_mercator +% Options for make_OGCM_SODA or make_OGCM_mercator %-------------------------------------------- % OGCM = 'SODA'; % Select OGCM: SODA or mercator @@ -430,6 +430,10 @@ ini_prefix = [CROCO_files_dir,'croco_ini_',OGCM,'_']; % generic initial file OGCM_prefix = [OGCM,'_']; % generic OGCM file name +% SODA +% ==== Description : 1/4 deg SODA v2.2.4 monthly means reanalysis (187101 -> 201012) + + % MERCATOR cases: % ============= mercator_type=1; % 1 --> 1/12 deg Mercator global reanalysis (1993 -> 202101) @@ -498,9 +502,9 @@ rmdepth = 2; itolap_a = 1; % before itolap_p = 1; % after % -%-------------------------- -% Options for make_bry_WKB -%-------------------------- +%---------------------------------------------------- +% Options for make_bry_WKB_ECMWF or make_bry_WKB_ERA5 +%---------------------------------------------------- % wkb_prefix=[CROCO_files_dir,'croco_wkb']; wkb_obc= [1 1 1 1]; @@ -510,6 +514,7 @@ wkb_obc= [1 1 1 1]; % 8 - Parameters for the forecast system % % --> select OGCM name above (mercator ...) +% --> select Aforc name (GFS/IFS) % --> don't forget to define in cppdefs.h: % - ROBUST_DIAG % - CLIMATOLOGY @@ -522,13 +527,8 @@ FRCST_dir = [FORC_DATA_DIR,'Forecast/']; % path to local OGCM data directory % % Number of hindcast / forecast days % -if strcmp(OGCM,'ECCO') - hdays=1; - fdays=6; -elseif strcmp(OGCM,'mercator') - hdays=1; - fdays=4; -end +hdays=1; +fdays=3; % % Local time= UTC + timezone %